Changeset 1103
 Timestamp:
 Oct 13, 2013 7:55:07 AM (9 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIrestrGUI.py
r1099 r1103 592 592 continue 593 593 else: 594 for res,name,id in Chains[chain][residue]: 595 try: 596 ipos = Atms.index(name) 597 ids[ipos] = id 598 except ValueError: 599 continue 600 for res,name,id in Chains[chain][residue+1]: 601 try: 602 ipos = pAtms.index(name) 603 ids[ipos] = id 604 except ValueError: 605 continue 594 try: 595 for res,name,id in Chains[chain][residue]: 596 try: 597 ipos = Atms.index(name) 598 ids[ipos] = id 599 except ValueError: 600 continue 601 for res,name,id in Chains[chain][residue+1]: 602 try: 603 ipos = pAtms.index(name) 604 ids[ipos] = id 605 except ValueError: 606 continue 607 except KeyError: 608 continue 606 609 if np.all(ids): 607 610 torsion = [list(ids),['1','1','1','1'],Name,Esd] … … 667 670 else: 668 671 for residue in residues[1:1]: 669 for res,name,id in Chains[chain][residue1]: 670 try: 671 ipos = mAtms.index(name) 672 ids[ipos] = id 673 except ValueError: 674 continue 675 for res,name,id in Chains[chain][residue+1]: 676 try: 677 ipos = pAtms.index(name) 678 ids[ipos] = id 679 except ValueError: 680 continue 681 for res,name,id in Chains[chain][residue]: 682 if Res == res: 672 try: 673 for res,name,id in Chains[chain][residue1]: 683 674 try: 684 ipos = Atms.index(name)675 ipos = mAtms.index(name) 685 676 ids[ipos] = id 686 677 except ValueError: 687 678 continue 688 if np.all(ids): 689 rama = [list(ids),['1','1','1','1','1'],Name,Esd] 690 if rama not in ramaRestData['Ramas']: 691 ramaRestData['Ramas'].append(rama) 692 ids = np.array([0,0,0,0,0]) 679 for res,name,id in Chains[chain][residue+1]: 680 try: 681 ipos = pAtms.index(name) 682 ids[ipos] = id 683 except ValueError: 684 continue 685 for res,name,id in Chains[chain][residue]: 686 if Res == res: 687 try: 688 ipos = Atms.index(name) 689 ids[ipos] = id 690 except ValueError: 691 continue 692 if np.all(ids): 693 rama = [list(ids),['1','1','1','1','1'],Name,Esd] 694 if rama not in ramaRestData['Ramas']: 695 ramaRestData['Ramas'].append(rama) 696 ids = np.array([0,0,0,0,0]) 697 except KeyError: 698 continue 693 699 macStr = macro.readline() 694 700 macro.close() … … 1813 1819 ramaRestData = restrData['Rama'] 1814 1820 UpdateRamaRestr(ramaRestData) 1815 G2plt.PlotRama(G2frame,phaseName,rama,ramaName)1821 wx.CallAfter(G2plt.PlotRama,G2frame,phaseName,rama,ramaName) 1816 1822 elif text == 'Chem. comp. restraints': 1817 1823 G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.RestraintMenu) 
trunk/GSASIIstrMath.py
r1101 r1103 591 591 refl[10] = atan2d(fbs[0],fas[0]) 592 592 593 def StructureFactor2(refList,G,hfx,pfx,SGData,calcControls,parmDict): 594 ''' Compute structure factors for all h,k,l for phase 595 puts the result, F^2, in each ref[8] in refList  want to do this for blocks of reflections 596 & not one by one. 597 input: 598 599 :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv] 600 :param np.array G: reciprocal metric tensor 601 :param str pfx: phase id string 602 :param dict SGData: space group info. dictionary output from SpcGroup 603 :param dict calcControls: 604 :param dict ParmDict: 605 606 ''' 607 twopi = 2.0*np.pi 608 twopisq = 2.0*np.pi**2 609 phfx = pfx.split(':')[0]+hfx 610 ast = np.sqrt(np.diag(G)) 611 Mast = twopisq*np.multiply.outer(ast,ast) 612 FFtables = calcControls['FFtables'] 613 BLtables = calcControls['BLtables'] 614 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 615 FF = np.zeros(len(Tdata)) 616 if 'N' in calcControls[hfx+'histType']: 617 FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam']) 618 else: 619 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 620 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 621 Uij = np.array(G2lat.U6toUij(Uijdata)) 622 bij = Mast*Uij.T 623 for refl in refList: 624 fbs = np.array([0,0]) 625 H = refl[:3] 626 SQ = 1./(2.*refl[4])**2 627 SQfactor = 4.0*SQ*twopisq 628 Bab = parmDict[phfx+'BabA']*np.exp(parmDict[phfx+'BabU']*SQfactor) 629 if not len(refl[1]): #no form factors 630 if 'N' in calcControls[hfx+'histType']: 631 refl[1] = G2el.getBLvalues(BLtables) 632 else: #'X' 633 refl[1] = G2el.getFFvalues(FFtables,SQ) 634 for i,El in enumerate(Tdata): 635 FF[i] = refl[1][El] 636 Uniq = refl[11] 637 phi = refl[12] 638 phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis]) 639 sinp = np.sin(phase) 640 cosp = np.cos(phase) 641 occ = Mdata*Fdata/len(Uniq) 642 biso = SQfactor*Uisodata 643 Tiso = np.where(biso<1.,np.exp(biso),1.0) 644 HbH = np.array([np.inner(h,np.inner(bij,h)) for h in Uniq]) 645 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 646 Tcorr = Tiso*Tuij 647 fa = np.array([(FF+FPBab)*occ*cosp*Tcorr,FPP*occ*sinp*Tcorr]) 648 fas = np.sum(np.sum(fa,axis=1),axis=1) #real 649 if not SGData['SGInv']: 650 fb = np.array([(FF+FPBab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr]) 651 fbs = np.sum(np.sum(fb,axis=1),axis=1) 652 fasq = fas**2 653 fbsq = fbs**2 #imaginary 654 refl[9] = np.sum(fasq)+np.sum(fbsq) 655 refl[10] = atan2d(fbs[0],fas[0]) 656 593 657 def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict): 594 658 'Needs a doc string' 659 twopi = 2.0*np.pi 660 twopisq = 2.0*np.pi**2 661 phfx = pfx.split(':')[0]+hfx 662 ast = np.sqrt(np.diag(G)) 663 Mast = twopisq*np.multiply.outer(ast,ast) 664 FFtables = calcControls['FFtables'] 665 BLtables = calcControls['BLtables'] 666 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 667 FF = np.zeros(len(Tdata)) 668 if 'N' in calcControls[hfx+'histType']: 669 FP = 0. 670 FPP = 0. 671 else: 672 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 673 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 674 Uij = np.array(G2lat.U6toUij(Uijdata)) 675 bij = Mast*Uij.T 676 dFdvDict = {} 677 dFdfr = np.zeros((len(refList),len(Mdata))) 678 dFdx = np.zeros((len(refList),len(Mdata),3)) 679 dFdui = np.zeros((len(refList),len(Mdata))) 680 dFdua = np.zeros((len(refList),len(Mdata),6)) 681 dFdbab = np.zeros((len(refList),2)) 682 for iref,refl in enumerate(refList): 683 H = np.array(refl[:3]) 684 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 685 SQfactor = 8.0*SQ*np.pi**2 686 dBabdA = np.exp(parmDict[phfx+'BabU']*SQfactor) 687 Bab = parmDict[phfx+'BabA']*dBabdA 688 for i,El in enumerate(Tdata): 689 FF[i] = refl[1][El] 690 Uniq = refl[11] 691 phi = refl[12] 692 phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:]) 693 sinp = np.sin(phase) 694 cosp = np.cos(phase) 695 occ = Mdata*Fdata/len(Uniq) 696 biso = SQfactor*Uisodata 697 Tiso = np.where(biso<1.,np.exp(biso),1.0) 698 HbH = np.inner(H,np.inner(bij,H)) 699 Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq]) 700 Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij]) 701 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 702 Tcorr = Tiso*Tuij 703 fot = (FF+FPBab)*occ*Tcorr 704 fotp = FPP*occ*Tcorr 705 fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp]) #non positions 706 fb = np.array([fot[:,np.newaxis]*sinp,fotp[:,np.newaxis]*sinp]) 707 708 fas = np.sum(np.sum(fa,axis=1),axis=1) 709 fbs = np.sum(np.sum(fb,axis=1),axis=1) 710 fax = np.array([fot[:,np.newaxis]*sinp,fotp[:,np.newaxis]*sinp]) #positions 711 fbx = np.array([fot[:,np.newaxis]*cosp,fot[:,np.newaxis]*cosp]) 712 #sum below is over Uniq 713 dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2) 714 dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2) 715 dfadui = np.sum(SQfactor*fa,axis=2) 716 dfadua = np.sum(Hij*fa[:,:,:,np.newaxis],axis=2) 717 dfadba = np.sum(cosp*(occ*Tcorr)[:,np.newaxis],axis=1) 718 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for 719 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq) 720 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1]) 721 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1]) 722 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1]) 723 dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 724 if not SGData['SGInv']: 725 dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom 726 dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2) 727 dfbdui = np.sum(SQfactor*fb,axis=2) 728 dfbdua = np.sum(Hij*fb[:,:,:,np.newaxis],axis=2) 729 dfbdba = np.sum(sinp*(occ*Tcorr)[:,np.newaxis],axis=1) 730 dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 731 dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 732 dFdui[iref] += 2.*(fbs[0]*dfbdui[0]fbs[1]*dfbdui[1]) 733 dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 734 dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 735 #loop over atoms  each dict entry is list of derivatives for all the reflections 736 for i in range(len(Mdata)): 737 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 738 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 739 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 740 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 741 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 742 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 743 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 744 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 745 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] 746 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] 747 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] 748 dFdvDict[pfx+'BabA'] = dFdbab.T[0] 749 dFdvDict[pfx+'BabU'] = dFdbab.T[1] 750 return dFdvDict 751 752 def StructureFactorDerv2(refList,G,hfx,pfx,SGData,calcControls,parmDict): 753 '''Needs a doc string  want to do this for blocks of reflections 754 & not one by one.''' 595 755 twopi = 2.0*np.pi 596 756 twopisq = 2.0*np.pi**2
Note: See TracChangeset
for help on using the changeset viewer.