Changeset 954


Ignore:
Timestamp:
Jun 18, 2013 3:53:19 PM (8 years ago)
Author:
vondreele
Message:

fix to Hessian LSQ to stop if lam > 10e5!
fix to rigid body quaternion derivatives

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r953 r954  
    130130                lam /= 10.
    131131                break
     132            if lam > 10.e5:
     133                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
     134                break
    132135        if (chisq0-chisq1)/chisq0 < ftol:
    133136            break
     
    350353        XYZ.append(atom[cx:cx+3])
    351354    return np.array(XYZ)
     355
     356def RotateRBXYZ(Bmat,Cart,oriQ):
     357    '''rotate & transform cartesian coordinates to crystallographic ones
     358    no translation applied. To be used for numerical derivatives
     359   
     360    :param type name: description
     361   
     362    :returns: type name: description
     363   
     364    '''
     365    ''' returns crystal coordinates for atoms described by RBObj
     366    '''
     367    XYZ = np.zeros_like(Cart)
     368    for i,xyz in enumerate(Cart):
     369        X = prodQVQ(oriQ,xyz)
     370        XYZ[i] = np.inner(Bmat,X)
     371    return XYZ
    352372
    353373def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
  • trunk/GSASIIstrMain.py

    r946 r954  
    4141   
    4242ateln2 = 8.0*math.log(2.0)
    43 DEBUG = False
     43DEBUG = True
    4444
    4545
  • trunk/GSASIIstrMath.py

    r953 r954  
    188188            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
    189189                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
    190             for iv in range(4):         #there is a problem with the Oa,Oi,Oj,Ok derivatives
     190            for iv in range(4):
    191191                Q[iv] -= dx
    192                 XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
     192                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
    193193                Q[iv] += 2.*dx
    194                 XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
     194                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
    195195                Q[iv] -= dx
    196196                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
     
    251251        for ia,atId in enumerate(RBObj['Ids']):
    252252            atNum = AtLookup[atId]
    253             dx = 0.0001
     253            dx = 0.00001
    254254            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
    255255                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
    256256            for iv in range(4):
    257257                Q[iv] -= dx
    258                 XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
     258                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
    259259                Q[iv] += 2.*dx
    260                 XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
     260                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
    261261                Q[iv] -= dx
    262262                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
Note: See TracChangeset for help on using the changeset viewer.