Changeset 3415


Ignore:
Timestamp:
May 31, 2018 3:47:44 PM (3 years ago)
Author:
vondreele
Message:

change testDeriv.py so it can be run in a bat file.
mag moment derivatives now close to correct & are workable

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r3414 r3415  
    860860        else:
    861861            fotp = FPP*Tcorr     
    862 #            GSASIIpath.IPyBreak()
    863862        if 'T' in calcControls[hfx+'histType']:
    864863            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     
    910909        dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \
    911910                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
    912 #        GSASIIpath.IPyBreak()
    913911        iBeg += blkSize
    914912#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
     
    982980    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
    983981    Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
    984     Kdata = np.inner(Gdata.T,uAmat).T*np.sqrt(nl.det(Ginv))/Mag     #Cartesian unit vectors
     982    VGi = np.sqrt(nl.det(Ginv))
     983    Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag     #Cartesian unit vectors
    985984    Uij = np.array(G2lat.U6toUij(Uijdata))
    986985    bij = Mast*Uij.T
     
    10751074    :returns: dict dFdvDict: dictionary of derivatives
    10761075    '''
     1076   
    10771077    g = nl.inv(G)
    10781078    ast = np.sqrt(np.diag(G))
     
    10971097    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
    10981098    dMdm = np.inner(Gdata.T,Ginv).T/Mag
     1099    Gones = np.ones_like(Gdata)
    10991100    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     1101    Gones = np.inner(Gones.T,SGMT).T
    11001102    if SGData['SGInv'] and not SGData['SGFixed']:
    11011103        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     1104        Gones = np.hstack((Gones,-Gones))       #inversion if any
    11021105    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
     1106    Gones = np.hstack([Gones for icen in range(Ncen)])        #dup over cell centering
    11031107    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
     1108    Gones = SGData['MagMom'][nxs,:,nxs]*Gones   #flip vectors according to spin flip
    11041109    Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
    11051110    VGi = np.sqrt(nl.det(Ginv))
    11061111    Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag       #make unit vectors in Cartesian space
    1107     dkdG = (np.inner(np.ones(3),uAmat)*VGi)[:,nxs,nxs]/Mag[nxs,:,:]
     1112#    Gones = np.ones_like(Gdata)
     1113    dkdG = (np.inner(Gones.T,uAmat).T*VGi)/Mag
    11081114    dkdm = dkdG-Kdata*dMdm[:,nxs,:]/Mag[nxs,:,:]
    11091115    dFdMx = np.zeros((nRef,mSize,3))
     
    11591165        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
    11601166        dqdk = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T     #Mxyz**2,Nref
    1161 #        NQ = np.where(np.abs(Q)>0.,1./np.abs(Q),0.)     #this sort of works esp for 1 axis moments
    11621167        NQ = np.sqrt(np.sum(Q*Q,axis=0))
    11631168        NQ = np.where(NQ > 0.,1./NQ,0.)
    11641169        dqdm = dqdk[:,:,:,nxs,nxs]*dkdm[:,nxs,nxs,:,:]   #Mxyz**2,Nref,Nops,Natms
    1165         dmx = NQ*Q*dMdm[:,nxs,nxs,:]
     1170        dmx = Q*dMdm[:,nxs,nxs,:]/2.
    11661171        dmx = dmx[nxs,:,:,:,:]+dqdm*Mag[nxs,nxs,nxs,:,:]
    11671172       
     
    11771182        dfadmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*cosm[nxs,nxs,:,:,:],axis=3)
    11781183        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
    1179         dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)    #deriv OK? not U12 & U23 in sarc
     1184        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK
    11801185        # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
    11811186        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
     
    11871192        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
    11881193        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
    1189         dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0)/len(SGT)          #problems
     1194        dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0)          #problems
    11901195        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
    11911196        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok
  • trunk/testDeriv.py

    r3369 r3415  
    2222
    2323import sys
     24import os
    2425import platform
    2526if '2' in platform.python_version_tuple()[0]:
     
    6970        self.testDerivPanel = wx.ScrolledWindow(self)       
    7071        self.plotNB = plot.PlotNotebook()
     72        self.testFile = ''
     73        arg = sys.argv
     74        if len(arg) > 1 and arg[1]:
     75            try:
     76                self.testFile = os.path.splitext(arg[1])[0]+u'.testDeriv'
     77            except:
     78                self.testFile = os.path.splitext(arg[1])[0]+'.testDeriv'
     79            self.TestRead()
     80            self.UpdateControls(None)
    7181       
    7282    def __init__(self, parent):
     
    94104            if dlg.ShowModal() == wx.ID_OK:
    95105                self.dirname = dlg.GetDirectory()
    96                 testFile = dlg.GetPath()
    97                 file = open(testFile,'rb')
    98                 if '2' in platform.python_version_tuple()[0]:
    99                     self.values = cPickle.load(file)
    100                     self.HistoPhases = cPickle.load(file)
    101                     (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
    102                     self.parmDict = cPickle.load(file)
    103                     self.varylist = cPickle.load(file)
    104                     self.calcControls = cPickle.load(file)
    105                     self.pawleyLookup = cPickle.load(file)
    106                 else:
    107                     self.values = cPickle.load(file,encoding='Latin-1')
    108                     self.HistoPhases = cPickle.load(file,encoding='Latin-1')
    109                     (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin-1')
    110                     self.parmDict = cPickle.load(file,encoding='Latin-1')
    111                     self.varylist = cPickle.load(file,encoding='Latin-1')
    112                     self.calcControls = cPickle.load(file,encoding='Latin-1')
    113                     self.pawleyLookup = cPickle.load(file,encoding='Latin-1')
    114                 self.use = [False for i in range(len(self.varylist+self.depVarList))]
    115                 self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.varylist+self.depVarList]
    116                 file.close()
    117                 groups,parmlist = G2mv.GroupConstraints(self.constrDict)
    118                 G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
     106                self.testFile = dlg.GetPath()
     107                self.TestRead()
    119108                self.UpdateControls(event)
    120                 print(G2mv.VarRemapShow(self.varylist))
    121                 print('Dependent Vary List:',self.depVarList)
    122109        finally:
    123110            dlg.Destroy()
     111           
     112    def TestRead(self):
     113        file = open(self.testFile,'rb')
     114        if '2' in platform.python_version_tuple()[0]:
     115            self.values = cPickle.load(file)
     116            self.HistoPhases = cPickle.load(file)
     117            (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
     118            self.parmDict = cPickle.load(file)
     119            self.varylist = cPickle.load(file)
     120            self.calcControls = cPickle.load(file)
     121            self.pawleyLookup = cPickle.load(file)
     122        else:
     123            self.values = cPickle.load(file,encoding='Latin-1')
     124            self.HistoPhases = cPickle.load(file,encoding='Latin-1')
     125            (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin-1')
     126            self.parmDict = cPickle.load(file,encoding='Latin-1')
     127            self.varylist = cPickle.load(file,encoding='Latin-1')
     128            self.calcControls = cPickle.load(file,encoding='Latin-1')
     129            self.pawleyLookup = cPickle.load(file,encoding='Latin-1')
     130        self.use = [False for i in range(len(self.varylist+self.depVarList))]
     131        self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.varylist+self.depVarList]
     132        file.close()
     133        groups,parmlist = G2mv.GroupConstraints(self.constrDict)
     134        G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
     135        print(G2mv.VarRemapShow(self.varylist))
     136        print('Dependent Vary List:',self.depVarList)
    124137           
    125138    def UpdateControls(self,event):
Note: See TracChangeset for help on using the changeset viewer.