Changeset 889 for trunk/GSASIIstruct.py

Timestamp:

Apr 19, 2013 3:38:07 PM (9 years ago)

Author:

vondreele

Message:

add torsion angle derivatives in RB
small fixes to other RB stuff

File:

• trunk/GSASIIstruct.py (modified) (12 diffs)

trunk/GSASIIstruct.py

@@ -814 +814 @@
         Pos = RBObj['Orig'][0]
         jrb = RRBIds.index(RBObj['RBId'])
+        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
         rbsx = str(irb)+':'+str(jrb)
         XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
+        for itors,tors in enumerate(RBObj['Torsions']):
+            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
+            orId,pvId = torData[itors][:2]
+            pivotVec = Cart[pvId]-Cart[orId]
+            QA = G2mth.AVdeg2Q(0.001,pivotVec)
+            QB = G2mth.AVdeg2Q(-0.001,pivotVec)
+            for ir in torData[itors][3]:
+                atNum = AtLookup[RBObj['Ids'][ir]]
+                rVec = Cart[ir]-Cart[pvId]
+                dRdT = np.inner(Bmat,(G2mth.prodQVQ(QA,rVec)-G2mth.prodQVQ(QB,rVec)))/.002
+                for ix in [0,1,2]:
+                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
         for ia,atId in enumerate(RBObj['Ids']):
             atNum = AtLookup[atId]
@@ -1287 +1300 @@
                     print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
                     coeffDict = itemRest['Coeff']
-                    for indx,ops,cofName,esd in enumerate(itemRest[rest]):
+                    for indx,ops,cofName,esd in itemRest[rest]:
                         AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
                         AtName = ''
@@ -1297 +1310 @@
                             tor = G2mth.getRestTorsion(XYZ,Amat)
                             restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
-                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(AtName[:-1],calc,obs,esd,(obs-calc)/esd,tor)
+                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
                         else:
                             phi,psi = G2mth.getRestRama(XYZ,Amat)
                             restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
-                            print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(AtName[:-1],calc,obs,esd,(obs-calc)/esd,phi,psi)
+                            print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])
                 elif name == 'ChemComp':
                     print >>pFile,'     atoms   mul*frac  factor     prod'
@@ -1469 +1482 @@
         valstr = '  values:'
         sigstr = '  esds  :'
+        if 'N' not in TLS:
+            print >>pFile,' Thermal motion:'
         if 'T' in TLS:
             for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
@@ -1513 +1528 @@
         if 'U' in TLS:
             name = pfx+rbfx+'U:'+rbsx
-            namstr = '  names :'+'%12s'%('U')
+            namstr = '  names :'+'%12s'%('Uiso')
             valstr = '  values:'+'%12.5f'%(parmDict[name])
             if name in sigDict:
@@ -1526 +1541 @@
         namstr = '  names :'
         valstr = '  values:'
-        sigstr = '  esds  :'        
+        sigstr = '  esds  :'
+        nTors = len(RBObj['Torsions'])
+        if nTors:
+            print >>pFile,' Torsions:'
+            for it in range(nTors):
+                name = pfx+'RBRTr;'+str(it)+':'+rbsx
+                namstr += '%12s'%('Tor'+str(it))
+                valstr += '%12.4f'%(parmDict[name])
+                if name in sigDict:
+                    sigstr += '%12.4f'%(sigDict[name])
+            print >>pFile,namstr
+            print >>pFile,valstr
+            print >>pFile,sigstr
                 
     def PrintSHtextureAndSig(textureData,SHtextureSig):
@@ -1626 +1653 @@
                 jrb = VRBIds.index(RBObj['RBId'])
                 rbsx = str(irb)+':'+str(jrb)
-                print >>pFile,' Vector rigid body '
+                print >>pFile,' Vector rigid body parameters:'
                 PrintRBObjPOAndSig('RBV',rbsx)
                 PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
@@ -1632 +1659 @@
                 jrb = RRBIds.index(RBObj['RBId'])
                 rbsx = str(irb)+':'+str(jrb)
-                print >>pFile,' Residue rigid body '
+                print >>pFile,' Residue rigid body parameters:'
                 PrintRBObjPOAndSig('RBR',rbsx)
                 PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
@@ -2683 +2710 @@
                 elif name in ['Torsion','Rama']:
                     coeffDict = itemRest['Coeff']
-                    for i,[indx,ops,cofName,esd] in enumerate(torsionList):
+                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
                         pNames.append(str(pId)+':'+name+':'+str(i))
                         XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
@@ -3981 +4008 @@
                         sig = ref[6]/(2.0*Fo)
                         if Fo/sig >= calcControls['minF/sig']:
-                            w = wtFactor/sig
+                            w = wtFactor/np.sqrt(sig)
                             for j,var in enumerate(varylist):
                                 if var in dFdvDict:
@@ -4064 +4091 @@
                         if ref[5]/ref[6] >= calcControls['minF/sig']:
                             w =  wtFactor/ref[6]
+                            print ref[:3],ref[5],ref[6],w
                             wdf[iref] = w*(ref[5]-ref[7])
                             for j,var in enumerate(varylist):
@@ -4077 +4105 @@
                             w = wtFactor/sig
                             wdf[iref] = w*(Fo-Fc)
+                            print ref[:3],ref[5],Fo,ref[6],sig,w
                             if Fo/sig >= calcControls['minF/sig']:
                                 for j,var in enumerate(varylist):