Changeset 3415
 Timestamp:
 May 31, 2018 3:47:44 PM (3 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIstrMath.py
r3414 r3415 860 860 else: 861 861 fotp = FPP*Tcorr 862 # GSASIIpath.IPyBreak()863 862 if 'T' in calcControls[hfx+'histType']: 864 863 fa = np.array([fot*cosp,np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) … … 910 909 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)])+ \ 911 910 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()913 911 iBeg += blkSize 914 912 # print 'derv time %.4f, nref %d, blkSize %d'%(time.time()time0,nRef,blkSize) … … 982 980 Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata #flip vectors according to spin flip * det(opM) 983 981 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 985 984 Uij = np.array(G2lat.U6toUij(Uijdata)) 986 985 bij = Mast*Uij.T … … 1075 1074 :returns: dict dFdvDict: dictionary of derivatives 1076 1075 ''' 1076 1077 1077 g = nl.inv(G) 1078 1078 ast = np.sqrt(np.diag(G)) … … 1097 1097 Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T]) 1098 1098 dMdm = np.inner(Gdata.T,Ginv).T/Mag 1099 Gones = np.ones_like(Gdata) 1099 1100 Gdata = np.inner(Gdata.T,SGMT).T #apply sym. ops. 1101 Gones = np.inner(Gones.T,SGMT).T 1100 1102 if SGData['SGInv'] and not SGData['SGFixed']: 1101 1103 Gdata = np.hstack((Gdata,Gdata)) #inversion if any 1104 Gones = np.hstack((Gones,Gones)) #inversion if any 1102 1105 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 1103 1107 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 1104 1109 Mag = np.tile(Mag[:,nxs],Nops).T #make Mag same length as Gdata 1105 1110 VGi = np.sqrt(nl.det(Ginv)) 1106 1111 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 1108 1114 dkdm = dkdGKdata*dMdm[:,nxs,:]/Mag[nxs,:,:] 1109 1115 dFdMx = np.zeros((nRef,mSize,3)) … … 1159 1165 Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK 1160 1166 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 moments1162 1167 NQ = np.sqrt(np.sum(Q*Q,axis=0)) 1163 1168 NQ = np.where(NQ > 0.,1./NQ,0.) 1164 1169 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. 1166 1171 dmx = dmx[nxs,:,:,:,:]+dqdm*Mag[nxs,nxs,nxs,:,:] 1167 1172 … … 1177 1182 dfadmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*cosm[nxs,nxs,:,:,:],axis=3) 1178 1183 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 sarc1184 dfadua = np.sum(Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2) #deriv OK 1180 1185 # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6) 1181 1186 dfbdfr = np.sum(fbm/occ,axis=2) #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem … … 1187 1192 dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok 1188 1193 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)#problems1194 dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0) #problems 1190 1195 dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0) #ok 1191 1196 dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0) #ok 
trunk/testDeriv.py
r3369 r3415 22 22 23 23 import sys 24 import os 24 25 import platform 25 26 if '2' in platform.python_version_tuple()[0]: … … 69 70 self.testDerivPanel = wx.ScrolledWindow(self) 70 71 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) 71 81 72 82 def __init__(self, parent): … … 94 104 if dlg.ShowModal() == wx.ID_OK: 95 105 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='Latin1') 108 self.HistoPhases = cPickle.load(file,encoding='Latin1') 109 (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin1') 110 self.parmDict = cPickle.load(file,encoding='Latin1') 111 self.varylist = cPickle.load(file,encoding='Latin1') 112 self.calcControls = cPickle.load(file,encoding='Latin1') 113 self.pawleyLookup = cPickle.load(file,encoding='Latin1') 114 self.use = [False for i in range(len(self.varylist+self.depVarList))] 115 self.delt = [max(abs(self.parmDict[name])*0.0001,1e6) 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() 119 108 self.UpdateControls(event) 120 print(G2mv.VarRemapShow(self.varylist))121 print('Dependent Vary List:',self.depVarList)122 109 finally: 123 110 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='Latin1') 124 self.HistoPhases = cPickle.load(file,encoding='Latin1') 125 (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin1') 126 self.parmDict = cPickle.load(file,encoding='Latin1') 127 self.varylist = cPickle.load(file,encoding='Latin1') 128 self.calcControls = cPickle.load(file,encoding='Latin1') 129 self.pawleyLookup = cPickle.load(file,encoding='Latin1') 130 self.use = [False for i in range(len(self.varylist+self.depVarList))] 131 self.delt = [max(abs(self.parmDict[name])*0.0001,1e6) 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) 124 137 125 138 def UpdateControls(self,event):
Note: See TracChangeset
for help on using the changeset viewer.