Changeset 3774


Ignore:
Timestamp:
Jan 3, 2019 9:32:48 AM (3 years ago)
Author:
vondreele
Message:

fix super indexing problem in transposeHKLF
fix reflection generation for incommensurate mag case in G2lattice & G2pwd
clean up non Fourier modulation calcs & remove analytic derivative stuff (now numeric)
fix uij derivative bug
work on incommensurate magnetic sturcture factors - not working yet
clean up of testDeriv - better choices for delt & force reflection regeneration before each derivative test

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r3766 r3774  
    74027402            dlg.Destroy()
    74037403        Super = data[1]['Super']
    7404         refList,badRefs = G2lat.transposeHKLF(Trans,Super,refList)
     7404        isup = 0
     7405        if Super:
     7406            isup = 1
     7407        refList,badRefs = G2lat.transposeHKLF(Trans,isup,refList)
    74057408        if len(badRefs):    #do I want to list badRefs?
    74067409            G2frame.ErrorDialog('Failed transformation','Matrix yields fractional hkl indices')
  • trunk/GSASIIlattice.py

    r3736 r3774  
    13441344def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
    13451345    'needs a doc string'
     1346    ifMag = False
     1347    if 'MagSpGrp' in SGData:
     1348        ifMag = True
    13461349    HKLs = []
    13471350    vec = np.array(Vec)
     
    13621365                if d >= dmin:
    13631366                    HKLM = np.array([h,k,l,dH])
    1364                     if G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData):
     1367                    if (G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData)) or ifMag:
    13651368                        HKLs.append([h,k,l,dH,d])   
    13661369    return HKLs
  • trunk/GSASIImath.py

    r3768 r3774  
    14371437    for iatm in range(Ax.shape[0]):
    14381438        nx = 0
    1439 #        if 'ZigZag' in waveTypes[iatm]:
    1440 #            nx = 1
    1441 #            Tmm = Ax[iatm][0][:2]                       
    1442 #            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
    1443 #            ZtauXt[iatm],ZtauXx[iatm] = posZigZagDerv(glTau,Tmm,XYZmax)
    1444 #        elif 'Block' in waveTypes[iatm]:
    1445 #            nx = 1
    1446 #            Tmm = Ax[iatm][0][:2]                       
    1447 #            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
    1448 #            ZtauXt[iatm],ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax)
     1439        if 'ZigZag' in waveTypes[iatm]:
     1440            nx = 1
     1441        elif 'Block' in waveTypes[iatm]:
     1442            nx = 1
    14491443        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
    14501444        if nx:   
    1451             StauX[iatm][:-nx] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
    1452             CtauX[iatm][:-nx] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
     1445            StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
     1446            CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
    14531447        else:
    14541448            StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
     
    15411535    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
    15421536    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS]
    1543    
    1544 def ModulationDerv2(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
    1545     '''
    1546     H: array refBlk x ops X hklt proj to hkl
    1547     HP: array refBlk x ops X hklt proj to hkl
    1548     Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
    1549     '''
    1550    
    1551     Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
    1552     dGdMfC = np.zeros(Mf)
    1553     dGdMfS = np.zeros(Mf)
    1554     Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
    1555     dGdMxC = np.zeros(Mx)
    1556     dGdMxS = np.zeros(Mx)
    1557     Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
    1558     dGdMuC = np.zeros(Mu)
    1559     dGdMuS = np.zeros(Mu)
    1560    
    1561     D = twopi*H[:,:,3,nxs]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
    1562     HdotX = twopi*np.inner(HP,Xmod)        #refBlk x ops x atoms X ngl
    1563     HdotXD = HdotX+D[:,:,nxs,:]
    1564     if nWaves[2]:
    1565         Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
    1566         HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
    1567         HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
    1568         HbH = np.exp(-np.sum(HuH,axis=-2)) #refBlk x ops x atoms x ngl; sum waves - OK vs Modulation version
    1569         part1 = -np.exp(-HuH)*Fmod    #refBlk x ops x atoms x waves x ngl
    1570         dUdAu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
    1571         dUdBu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
    1572         dGdMuCa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
    1573         dGdMuCb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
    1574         dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
    1575         dGdMuSa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
    1576         dGdMuSb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
    1577         dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
    1578     else:
    1579         HbH = np.ones_like(HdotX)
    1580     dHdXA = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
    1581     dHdXB = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,nxs,:,:,:,:]    #ditto - cos waves
    1582 # ops x atoms x waves x 2xyz - real part - good
    1583     dGdMxCa = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
    1584     dGdMxCb = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
    1585     dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
    1586 # ops x atoms x waves x 2xyz - imag part - good
    1587     dGdMxSa = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
    1588     dGdMxSb = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
    1589     dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
    1590 # ZigZag/Block waves - problems here?
    1591     dHdXZt = -twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
    1592     dHdXZx = twopi*HP[:,:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
    1593     dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
    1594     dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    1595     dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
    1596     dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
    1597     dGdMzSx = np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    1598     dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
    1599 #    GSASIIpath.IPyBreak()
    1600     return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
    1601    
     1537       
    16021538def posFourier(tau,psin,pcos):
    16031539    A = np.array([ps[:,nxs]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
  • trunk/GSASIIphsGUI.py

    r3770 r3774  
    13421342                            atom[-1]['SS1'][parm] = [wType,]+list(atom[-1]['SS1'][parm])
    13431343                    del atom[-1]['SS1']['waveType']
     1344        else:
     1345            generalData['Super'] = 0
    13441346        if 'Modulated' not in generalData:
    13451347            generalData['Modulated'] = False
     
    15971599                            generalData['SSGData'] = G2spc.SSpcGroup(generalData['SGData'],generalData['SuperSg'])[1]
    15981600                            if 'SuperVec' not in generalData:
    1599                                 generalData['Super'] = 1
     1601                                generalData['Super'] = True
    16001602                                generalData['SuperVec'] = [[0.,0.,0.],False,4]
    16011603                                generalData['SSGData'] = {}
     
    50895091                            if np.any(CSI[0][ival]):
    50905092                                minmax = [-0.2,0.2]
     5093                                if Stype == 'Smag':
     5094                                    minmax = [-20.,20.]
    50915095                                if waveTyp in ['ZigZag','Block','Crenel'] and not iwave and ival < 2:
    50925096                                    if not ival:
  • trunk/GSASIIpwd.py

    r3712 r3774  
    10641064    SSdH = [vec*h for h in range(-maxH,maxH+1)]
    10651065    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
     1066    ifMag = False
     1067    if 'MagSpGrp' in SGData:
     1068        ifMag = True
    10661069    for h,k,l,d in HKL:
    10671070        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
     
    10751078                if d >= dmin:
    10761079                    HKLM = np.array([h,k,l,dH])
    1077                     if G2spc.checkSSextc(HKLM,SSGData):
     1080                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
    10781081                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
    10791082    return G2lat.sortHKLd(HKLs,True,True,True)
  • trunk/GSASIIstrMath.py

    r3773 r3774  
    850850        biso = -SQfactor*Uisodata[:,nxs]
    851851        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
    852         HbH = np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    853         Tuij = np.where(HbH<1.,np.exp(-HbH),1.0).T
     852        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     853        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    854854        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
    855855        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
     
    14971497            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
    14981498            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
    1499         Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
     1499        Bab = 0.
     1500        if phfx+'BabA' in parmDict:
     1501            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
    15001502        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    15011503        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
     
    15081510        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    15091511        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
     1512        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x refBlk x sym X atoms
     1513        fams = np.zeros(32)
     1514        fbms = np.zeros(32)
    15101515
    15111516        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: mag math here??
     
    15161521            else:
    15171522                mphase = phase
    1518             mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
     1523            mphase = np.array([mphase+twopi*np.inner(cen,HP.T)[:,nxs,nxs] for cen in SGData['SGCen']])
    15191524            mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
    15201525            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    15211526            cosm = np.cos(mphase)                               #ditto
    1522             HM = np.inner(Bmat,H)                             #put into cartesian space
     1527            HM = np.inner(Bmat,HP.T)                             #put into cartesian space
    15231528            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    15241529            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
     
    15261531            fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    15271532            fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    1528             fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
    1529             fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    1530             refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    1531             refl.T[7] = np.copy(refl.T[9])               
    1532             refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     1533            fagm = fam*GfpuA[0]-fbm*GfpuA[1]   #real; 2 x refBlk x sym x atoms
     1534            fbgm = fbm*GfpuA[0]+fam*GfpuA[1]
     1535            fams = np.sum(np.sum(fagm,axis=-1),axis=-1)                          #xyz,Nref
     1536            fbms = np.sum(np.sum(fbgm,axis=-1),axis=-1)                          #ditto
    15331537
    15341538
     1539        if 'T' in calcControls[hfx+'histType']:
     1540            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     1541            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    15351542        else:
    1536             if 'T' in calcControls[hfx+'histType']:
    1537                 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    1538                 fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1539             else:
    1540                 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    1541                 fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1542             GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x refBlk x sym X atoms
    1543             fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
    1544             fbg = fb*GfpuA[0]+fa*GfpuA[1]
    1545             fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
    1546             fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
    1547             if 'P' in calcControls[hfx+'histType']:
    1548                 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1549                 refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    1550             else:
    1551                 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1552                 refl.T[8] = np.copy(refl.T[10])               
    1553                 refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1543            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     1544            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1545        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
     1546        fbg = fb*GfpuA[0]+fa*GfpuA[1]
     1547        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
     1548        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
     1549        refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2+np.sum(fams,axis=0)**2+np.sum(fbms,axis=0)**2    #square of sums
     1550        refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1551        if 'P' not in calcControls[hfx+'histType']:
     1552            refl.T[8] = np.copy(refl.T[10])               
    15541553        iBeg += blkSize
    15551554#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
     
    17731772        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
    17741773        SQfactor = 8.0*SQ*np.pi**2
    1775         dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
    1776         Bab = parmDict[phfx+'BabA']*dBabdA
     1774        Bab = 0.0
     1775        if phfx+'BabA' in parmDict:
     1776            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1777            Bab = parmDict[phfx+'BabA']*dBabdA
    17771778        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    17781779        FF = refDict['FF']['FF'][iref].T[Tindx]
     
    18641865            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
    18651866            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
    1866         dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    1867             2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1867        if phfx+'BabA' in parmDict:
     1868            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1869                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    18681870        #loop over atoms - each dict entry is list of derivatives for all the reflections
    18691871        if not iref%100 :
  • trunk/testDeriv.py

    r3711 r3774  
    4444    return testDeriv(parent)
    4545   
    46 [wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS,
    47 ] = [wx.NewId() for _init_coll_File_Items in range(3)]
     46[wxID_FILEEXIT, wxID_FILEOPEN, wxID_MAKEPLOTS, wxID_CLEARSEL,
     47] = [wx.NewId() for _init_coll_File_Items in range(4)]
    4848
    4949def FileDlgFixExt(dlg,file):            #this is needed to fix a problem in linux wx.FileDialog
     
    6262        self.File.Append(wxID_FILEOPEN,'Open testDeriv file','Open testDeriv')
    6363        self.File.Append(wxID_MAKEPLOTS,'Make plots','Make derivative plots')
     64        self.File.Append(wxID_CLEARSEL,'Clear selections')
    6465        self.File.Append(wxID_FILEEXIT,'Exit','Exit from testDeriv')
    65         self.Bind(wx.EVT_MENU, self.OnTestRead, id=wxID_FILEOPEN)
     66        self.Bind(wx.EVT_MENU,self.OnTestRead, id=wxID_FILEOPEN)
    6667        self.Bind(wx.EVT_MENU,self.OnMakePlots,id=wxID_MAKEPLOTS)
    67         self.Bind(wx.EVT_MENU, self.OnFileExit, id=wxID_FILEEXIT)
     68        self.Bind(wx.EVT_MENU,self.ClearSelect,id=wxID_CLEARSEL)
     69        self.Bind(wx.EVT_MENU,self.OnFileExit, id=wxID_FILEEXIT)
    6870        self.testDerivMenu.Append(menu=self.File, title='File')
    6971        self.SetMenuBar(self.testDerivMenu)
     
    9597            self.dataFrame.Destroy()
    9698        self.Close()
     99       
     100    def ClearSelect(self,event):
     101        self.use = [False for i in range(len(self.names))]
     102        self.UpdateControls(event)
    97103
    98104    def OnTestRead(self,event):
     
    128134            self.calcControls = cPickle.load(file,encoding='Latin-1')
    129135            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]
     136        self.names = self.varylist+self.depVarList
     137        self.use = [False for i in range(len(self.names))]
     138        self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.names]
     139        for iname,name in enumerate(self.names):
     140            if name.split(':')[-1] in ['Shift','DisplaceX','DisplaceY',]:
     141                self.delt[iname] = 0.1
    132142        file.close()
    133143        msg = G2mv.EvaluateMultipliers(self.constrDict,self.parmDict)
     
    159169        self.testDerivPanel.DestroyChildren()
    160170        ObjInd = {}
    161         varylist = self.varylist
    162         depVarList = self.depVarList
     171        names = self.names
    163172        use = self.use
    164173        delt = self.delt
    165174        mainSizer = wx.FlexGridSizer(0,8,5,5)
    166         for id,[ck,name,d] in enumerate(zip(use,varylist+depVarList,delt)):
     175        for id,[ck,name,d] in enumerate(zip(use,names,delt)):
    167176            useVal = wx.CheckBox(self.testDerivPanel,label=name)
    168177            useVal.SetValue(ck)
     
    202211        def test2(name,delt,doProfile):
    203212            Title = 'derivatives test for '+name
    204             varyList = self.varylist+self.depVarList
     213            names = self.names
    205214            hplot = self.plotNB.add(Title).gca()
    206215            if doProfile:
    207216                pr = cProfile.Profile()
    208217                pr.enable()
     218            #regenerate minimization fxn
     219            G2stMth.errRefine(self.values,self.HistoPhases,
     220                self.parmDict,self.varylist,self.calcControls,
     221                self.pawleyLookup,None)
    209222            dMdV = G2stMth.dervRefine(self.values,self.HistoPhases,self.parmDict,
    210                 varyList,self.calcControls,self.pawleyLookup,None)
     223                names,self.calcControls,self.pawleyLookup,None)
    211224            if doProfile:
    212225                pr.disable()
     
    217230                print('Profiler of '+name+' derivative calculation; top 50% of routines:')
    218231                print(s.getvalue())
    219             M2 = dMdV[varyList.index(name)]
     232            M2 = dMdV[names.index(name)]
    220233            hplot.plot(M2,'b',label='analytic deriv')
    221             mmin = np.min(dMdV[varyList.index(name)])
    222             mmax = np.max(dMdV[varyList.index(name)])
     234            mmin = np.min(dMdV[names.index(name)])
     235            mmax = np.max(dMdV[names.index(name)])
    223236            print('parameter:',name,self.parmDict[name],delt,mmin,mmax)
    224237            if name in self.varylist:
    225238                self.values[self.varylist.index(name)] -= delt
    226239                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
    227                     varyList,self.calcControls,self.pawleyLookup,None)
     240                    names,self.calcControls,self.pawleyLookup,None)
    228241                self.values[self.varylist.index(name)] += 2.*delt
    229242                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
    230                     varyList,self.calcControls,self.pawleyLookup,None)
     243                    names,self.calcControls,self.pawleyLookup,None)
    231244                self.values[self.varylist.index(name)] -= delt
    232245            elif name in self.depVarList:   #in depVarList
     
    236249                self.parmDict[name] -= delt
    237250                M0 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
    238                     varyList,self.calcControls,self.pawleyLookup,None)
     251                    names,self.calcControls,self.pawleyLookup,None)
    239252                self.parmDict[name] += 2.*delt
    240253                M1 = G2stMth.errRefine(self.values,self.HistoPhases,self.parmDict,
    241                     varyList,self.calcControls,self.pawleyLookup,None)
     254                    names,self.calcControls,self.pawleyLookup,None)
    242255                self.parmDict[name] -= delt   
    243256            Mn = (M1-M0)/(2.*abs(delt))
Note: See TracChangeset for help on using the changeset viewer.