Changeset 820


Ignore:
Timestamp:
Jan 3, 2013 4:03:22 PM (10 years ago)
Author:
vondreele
Message:

restraints development - almost done
works for bonds, angles, planes & chiral volumes
not torsions or Rama.. yet

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/Exercises/sucrose/SUCROSE.EXP

    r534 r820  
    5353    HSTRY 51 EXPEDT  Win32  Sep 30 16:25:23 2007             L   O             
    5454    HSTRY 52 EXPEDT  Win32  Oct 19 09:36:53 2011             LA                 
     55    HSTRY 53 EXPEDT  Win32  Jan 03 10:20:47 2013 P  H        LA LO             
     56    HSTRY 54 POWPREF Win32  Jan 03 10:20:58 2013                               
     57    HSTRY 55 EXPEDT  Win32  Jan 03 12:11:16 2013             L S  D             
     58    HSTRY 56 EXPEDT  Win32  Jan 03 10:23:42 2013             L   O             
    5559  DSGL CDAT1  DRAD ARAD NOFO                                                   
    5660  DSGL CDAT2  DRAD ARAD NOFO                                                   
    … …  
    7680 AFAC  O_XF1  0.093 0.072 0.049 0.011 0.006 0.121 0.063-0.002-0.003-0.003       
    7781 AFAC  O_XF2  0.073 0.052 0.032 0.006 0.004 0.106 0.044 0.000 0.000 0.000       
    78  EXPR  HTYP1  PXC                                                              
     82 EXPR  HTYP1  PXC  RSN                                                         
    7983 EXPR  NATYP    3                                                               
    80  EXPR  NHST     1                                                               
     84 EXPR  NHST     2                                                               
    8185 EXPR ATYP 1  O          11      1.09      0.89 1.45   96 1.40                 
    8286 EXPR ATYP 2  C          12      1.12      0.92 1.65    3                       
    … …  
    8892CRS1    PNAM  SUCROSE                                                           
    8993CRS1   NATOM   45                                                               
    90 CRS1  ABC    10.863168  8.703946  7.757701    Y    0                           
     94CRS1  ABC    10.863168  8.703946  7.757701    N    0                           
    9195CRS1  ABCSIG  0.000153  0.000120  0.000107                                     
    92 CRS1  ANGLES   90.0000  102.9426   90.0000                                     
     96CRS1  ANGLES  90.00000 102.94260  90.00000                                     
    9397CRS1  ANGSIG    0.0000    0.0005    0.0000                                     
    9498CRS1  AT  1A  O        -0.164254 -0.157852  0.607290  1.000000O1         2 000 
    95 CRS1  AT  1B  0.040561                                                    I XU 
     99CRS1  AT  1B  0.040561                                                    I    
    96100CRS1  AT  2A  O        -0.231485 -0.061788  0.265443  1.000000O2         2 000 
    97 CRS1  AT  2B  0.057560                                                    I XU 
     101CRS1  AT  2B  0.057560                                                    I    
    98102CRS1  AT  3A  O        -0.295166  0.245759  0.317444  1.000000O3         2 000 
    99 CRS1  AT  3B  0.109726                                                    I XU 
     103CRS1  AT  3B  0.109726                                                    I    
    100104CRS1  AT  4A  O        -0.340114  0.326551  0.654647  1.000000O4         2 000 
    101 CRS1  AT  4B  0.138454                                                    I XU 
     105CRS1  AT  4B  0.138454                                                    I    
    102106CRS1  AT  5A  O        -0.392465 -0.103270  0.623136  1.000000O5         2 000 
    103 CRS1  AT  5B  0.028472                                                    I XU 
     107CRS1  AT  5B  0.028472                                                    I    
    104108CRS1  AT  6A  O        -0.576481  0.051765  0.726621  1.000000O6         2 000 
    105 CRS1  AT  6B  0.047865                                                    I XU 
     109CRS1  AT  6B  0.047865                                                    I    
    106110CRS1  AT  7A  O        -0.026064 -0.269992  0.384540  1.000000O11        2 000 
    107 CRS1  AT  7B  0.067787                                                    I XU 
     111CRS1  AT  7B  0.067787                                                    I    
    108112CRS1  AT  8A  O        -0.219269 -0.443970  0.672126  1.000000O21        2 000 
    109 CRS1  AT  8B  0.053494                                                    I XU 
     113CRS1  AT  8B  0.053494                                                    I    
    110114CRS1  AT  9A  O         0.067388 -0.191789  0.783754  1.000000O31        2 000 
    111 CRS1  AT  9B  0.061425                                                    I XU 
     115CRS1  AT  9B  0.061425                                                    I    
    112116CRS1  AT 10A  O         0.018320 -0.397973  1.080663  1.000000O41        2 000 
    113 CRS1  AT 10B  0.111212                                                    I XU 
     117CRS1  AT 10B  0.111212                                                    I    
    114118CRS1  AT 11A  O        -0.328317 -0.287414  0.964056  1.000000O61        2 000 
    115 CRS1  AT 11B  0.096466                                                    I XU 
     119CRS1  AT 11B  0.096466                                                    I    
    116120CRS1  AT 12A  C        -0.299305 -0.166713  0.523976  1.000000C1         2 000 
    117 CRS1  AT 12B  0.021275                                                    I XU 
     121CRS1  AT 12B  0.021275                                                    I    
    118122CRS1  AT 13A  C        -0.333469 -0.017226  0.345055  1.000000C2         2 000 
    119 CRS1  AT 13B  0.019679                                                    I XU 
     123CRS1  AT 13B  0.019679                                                    I    
    120124CRS1  AT 14A  C        -0.267194  0.121736  0.461451  1.000000C3         2 000 
    121 CRS1  AT 14B  0.018025                                                    I XU 
     125CRS1  AT 14B  0.018025                                                    I    
    122126CRS1  AT 15A  C        -0.374575  0.172912  0.565194  1.000000C4         2 000 
    123 CRS1  AT 15B  0.000268                                                    I XU 
     127CRS1  AT 15B  0.000268                                                    I    
    124128CRS1  AT 16A  C        -0.349910  0.106966  0.716598  1.000000C5         2 000 
    125 CRS1  AT 16B  0.032578                                                    I XU 
     129CRS1  AT 16B  0.032578                                                    I    
    126130CRS1  AT 17A  C        -0.446764  0.038868  0.802935  1.000000C6         2 000 
    127 CRS1  AT 17B  0.022298                                                    I XU 
     131CRS1  AT 17B  0.022298                                                    I    
    128132CRS1  AT 18A  C        -0.102176 -0.408149  0.454657  1.000000C11        2 000 
    129 CRS1  AT 18B  0.002823                                                    I XU 
     133CRS1  AT 18B  0.002823                                                    I    
    130134CRS1  AT 19A  C        -0.117978 -0.322890  0.643846  1.000000C21        2 000 
    131 CRS1  AT 19B -0.019107                                                    I XU 
     135CRS1  AT 19B -0.019107                                                    I    
    132136CRS1  AT 20A  C        -0.014363 -0.345739  0.794304  1.000000C31        2 000 
    133 CRS1  AT 20B  0.095309                                                    I XU 
     137CRS1  AT 20B  0.095309                                                    I    
    134138CRS1  AT 21A  C        -0.076372 -0.348567  0.921647  1.000000C41        2 000 
    135 CRS1  AT 21B  0.001039                                                    I XU 
     139CRS1  AT 21B  0.001039                                                    I    
    136140CRS1  AT 22A  C        -0.197330 -0.436359  0.845452  1.000000C51        2 000 
    137 CRS1  AT 22B  0.047428                                                    I XU 
     141CRS1  AT 22B  0.047428                                                    I    
    138142CRS1  AT 23A  C        -0.285790 -0.475785  0.972198  1.000000C61        2 000 
    139 CRS1  AT 23B  0.002722                                                    I XU 
     143CRS1  AT 23B  0.002722                                                    I    
    140144CRS1  AT 24A  H        -0.334700 -0.254900  0.461200  1.000000H1         2 000 
    141145CRS1  AT 24B  0.029000                                                    I     
    … …  
    192196HAP1 1NAXIS     1                                                               
    193197HAP1 1PHSFR      1.0000        N    0                                           
    194 HAP1 1PRCF      5   17     0.005      YY YY  Y YY                               
     198HAP1 1PRCF      5   17   0.00500    0NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN       
    195199HAP1 1PRCF 1   0.000000E+00  -0.998165E+00   0.202298E+01   0.000000E+00       
    196200HAP1 1PRCF 2   0.985688E+00   0.209294E+01   0.100000E-02   0.100000E-02       
    … …  
    200204HAP1 1PREFO1  1.000000  1.000000  0.000000  0.000000  1.000000   NN    0    0   
    201205HAP1 1RADDAM    0.00000E+00    N    0      0.00                                 
    202 HST  1  HFIL  C:\RBVD\sucrose\sucrose_1sec_A_0.GSA                              
    203 HST  1  HNAM  F:\Oct_2003\disk1\sucrose_std_745mm_0.mar3450: Azimuth/2-theta   
    204 HST  1  IFIL  C:\RBVD\sucrose\inst_xry.prm                                     
     206HST  1  HFIL  D:\pyGSAS\Exercises\sucrose\11bmb_8716.fxye                       
     207HST  1  HNAM  /data/jul12/11bmb_8716.mda JUL 26, 2012 15:09:18.910830          
     208HST  1  IFIL  D:\pyGSAS\Exercises\sucrose\11bmb_8716.prm                       
    205209HST  1 BANK     1                                                               
    206210HST  1 BIGFO   0.166072E+06                                                     
    207 HST  1 CHANS      1386         1      3613  30744340      5000                 
     211HST  1 CHANS      1501         1     22405  87900360     49493    0    0       
    208212HST  1 EPHAS    0    0    0    0    0    0    0    0    0    Y                 
    209 HST  1 ICONR 0.6194280 0.0000000  -0.13099    Z    0   0.95000    0   0.50000   
    210 HST  1 ICONS 0.6194280 0.0000000  -1.84763    Z    0   0.95000    0   0.50000   
     213HST  1 ICONR 0.4133190    0.0000    0.0000               0.990    0     0.      
     214HST  1 ICONS 0.4133190    0.0000    0.0000               0.990    0     0.      
    211215HST  1 IRAD     0                                                               
    212 HST  1 MAXRF    8                                                               
     216HST  1 MAXRF   62                                                               
    213217HST  1 NEXC     2                                                               
     218HST  1 NFOBS  855    0    0    0    0    0    0    0    0                       
    214219HST  1 NPHAS    1    0    0    0    0    0    0    0    0                       
    215 HST  1 NREF   133    1.8765    Y    Y                                           
     220HST  1 NREF   888    0.9777                                                     
    216221HST  1 R-FAC  123   0.01160   5.233365E+04                                     
    217222HST  1 RPOWD    0.0234    0.0196   3613  24879.                                 
    218 HST  1 TRNGE   7.00092  18.91080                                               
    219 HST  1 WREXP    0.0089                                                         
     223HST  1 TRNGE   2.00100  24.40500                                               
     224HST  1 WREXP    0.0354                                                         
    220225HST  1ABSCOR   0.000000E+00   0.000000E+00    N    0                           
    221226HST  1BAKGD     7   20    Y    0    Y                                           
    … …  
    226231HST  1BAKGD5   0.112117E+04   0.104554E+04   0.102951E+04   0.102649E+04       
    227232HST  1CHI         0.00                                                         
    228 HST  1EXC  1     0.000     7.000                                               
    229 HST  1EXC  2    19.000  1000.000                                               
     233HST  1EXC  1     0.000     2.000                                               
     234HST  1EXC  2    24.406  1000.000                                               
    230235HST  1EXMNMX   1.00000   1.00000                                               
    231236HST  1FFANS1  O            0.000     0.000   NN    0                           
    … …  
    233238HST  1FFANS3  H            0.000     0.000   NN    0                           
    234239HST  1HSCALE     54.989        Y    0                                           
    235 HST  1I HEAD  DUMMY INCIDENT SPECTRUM FOR X-RAY DIFFRACTOMETER                 
     240HST  1I HEAD  11-BM profile SRM 660a (LaB6) fit (Feb 2009)                     
    236241HST  1I ITYP    0    0.0000  180.0000         1                                 
     242HST  1INAME   APS11BM                                                           
    237243HST  1MNREF     0    1.9792                                                     
    238244HST  1ODMNMX   1.00000   1.00000                                               
    239 HST  1PRCF1     5    8     0.005                                               
    240 HST  1PRCF11   0.000000E+00   0.418369E+01   0.162835E+01   0.000000E+00       
    241 HST  1PRCF12   0.000000E+00   0.000000E+00   0.010000E-01   0.010000E-01       
     245HST  1PRCF1     3   19   0.00100                                               
     246HST  1PRCF11       1.163000      -0.126000       0.063000   0.000000E+00       
     247HST  1PRCF12       0.173000       0.000000       0.001700       0.001700       
     248HST  1PRCF13   0.000000E+00   0.000000E+00   0.000000E+00   0.000000E+00       
     249HST  1PRCF14   0.000000E+00   0.000000E+00   0.000000E+00   0.000000E+00       
     250HST  1PRCF15   0.000000E+00   0.000000E+00   0.000000E+00                       
    242251HST  1RSP    0.0308 0.0199 0.0274 0.0097 0.0113 0.0094 0.0119 0.0086 0.0106     
    243252HST  1RSPA   0.0455 0.0310 0.0354 0.0146 0.0124 0.0091 0.0123 0.0092 0.0102     
    … …  
    247256HST  1TRMNMX   1.00000   1.00000                                               
    248257ZZZZZZZZZZZZ  Last EXP file record                                             
     258    HSTRY 57 POWPREF Win32  Jan 03 10:23:52 2013                               
     259    HSTRY 58 GENLES  Win32  Jan 03 10:24:02 2013 Sdsq= 0.580E+07 S/E= 0.260E-03
     260    HSTRY 59 EXPEDT  Win32  Jan 03 10:25:59 2013             L   O             
     261    HSTRY 60 POWPREF Win32  Jan 03 10:26:05 2013                               
     262    HSTRY 61 GENLES  Win32  Jan 03 10:26:11 2013 Sdsq= 0.338E+07 S/E= 0.306E-04
     263    HSTRY 62 EXPEDT  Win32  Jan 03 10:27:27 2013             L   O             
     264    HSTRY 63 GENLES  Win32  Jan 03 10:27:34 2013 Sdsq= 0.108E+07 S/E=  55.2     
     265    HSTRY 64 EXPEDT  Win32  Jan 03 10:28:36 2013             LA                 
     266    HSTRY 65 POWPREF Win32  Jan 03 10:28:39 2013                               
     267    HSTRY 66 GENLES  Win32  Jan 03 10:28:52 2013 Sdsq= 0.544E+06 S/E=  4.70     
     268    HSTRY 67 EXPEDT  Win32  Jan 03 10:30:47 2013             LA  O             
     269    HSTRY 68 GENLES  Win32  Jan 03 10:31:00 2013 Sdsq= 0.256E+06 S/E=  110.     
     270    HSTRY 69 EXPEDT  Win32  Jan 03 10:35:34 2013             LAS  D             
     271HST  2  HNAM  Bond distance restraints                                         
     272HST  2 FACTR    1.00000                                                         
     273HST  2 NBNDS   22                                                               
     274HST  2BD0001  1   12   24  1  0  0  0  0 1.090 0.010                           
     275HST  2BD0002  1   13   25  1  0  0  0  0 1.090 0.010                           
     276HST  2BD0003  1   14   26  1  0  0  0  0 1.090 0.010                           
     277HST  2BD0004  1   15   27  1  0  0  0  0 1.090 0.010                           
     278HST  2BD0005  1   16   28  1  0  0  0  0 1.090 0.010                           
     279HST  2BD0006  1   17   29  1  0  0  0  0 1.090 0.010                           
     280HST  2BD0007  1   17   38  1  0  0  0  0 1.090 0.010                           
     281HST  2BD0008  1   18   31  1  0  0  0  0 1.090 0.010                           
     282HST  2BD0009  1   20   33  1  0  0  0  0 1.090 0.010                           
     283HST  2BD000A  1   21   35  1  0  0  0  0 1.090 0.010                           
     284HST  2BD000B  1   22   37  1  0  0  0  0 1.090 0.010                           
     285HST  2BD000C  1   23   39  1  0  0  0  0 1.090 0.010                           
     286HST  2BD000D  1   23   45  1  0  0  0  0 1.090 0.010                           
     287HST  2BD000E  1    2   32  1  0  0  0  0 1.090 0.010                           
     288HST  2BD000F  1    3   34  1  0  0  0  0 1.090 0.010                           
     289HST  2BD0010  1    4   36  1  0  0  0  0 1.090 0.010                           
     290HST  2BD0011  1    6   40  1  0  0  0  0 1.090 0.010                           
     291HST  2BD0012  1    7   41  1  0  0  0  0 1.090 0.010                           
     292HST  2BD0013  1    9   42  1  0  0  0  0 1.090 0.010                           
     293HST  2BD0014  1   10   43  1  0  0  0  0 1.090 0.010                           
     294HST  2BD0015  1   11   44  1  0  0  0  0 1.090 0.010                           
     295    HSTRY 70 GENLES  Win32  Jan 03 10:35:53 2013 Sdsq= 0.209E+06 S/E=  129.     
     296    HSTRY 71 EXPEDT  Win32  Jan 03 10:40:04 2013             L S  D             
     297 REFN RESTR      22  27.197                                                     
     298HST  2 RMSD    22    0.0004    27.1966                                         
     299HST  2BD0016  1   18   30  1  0  0  0  0 1.090 0.010                           
     300    HSTRY 72 POWPREF Win32  Jan 03 10:40:08 2013                               
     301    HSTRY 73 GENLES  Win32  Jan 03 10:40:27 2013 Sdsq= 0.204E+06 S/E=  90.4     
     302    HSTRY 74 EXPEDT  Win32  Jan 03 11:35:16 2013             L S  D             
     303    HSTRY 75 GENLES  Win32  Jan 03 11:35:35 2013 Sdsq= 0.208E+06 S/E=  12.9     
     304    HSTRY 76 EXPEDT  Win32  Jan 03 11:36:50 2013             L S  D             

  • trunk/GSASIImath.py

    r818 r820  
    186186            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
    187187    return Items
     188   
     189def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
     190    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
     191    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
     192    XYZ = []
     193    for ind in indx:
     194        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
     195        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
     196        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
     197    return XYZ
    188198       
    189199def getMass(generalData):
    … …  
    223233    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
    224234   
    225 def getRestDeriv(Func,XYZ,Amat):
    226     deriv = np.array((3,len(XYZ)))
     235def getRestDeriv(Func,XYZ,Amat,ops,SGData):
     236    deriv = np.zeros((len(XYZ),3))
    227237    dx = 0.00001
    228238    for j,xyz in enumerate(XYZ):
    229         for i,x in enumerate(xyz):
    230             x += dx
    231             d1 = Func(XYZ,Amat)
    232             x -= 2*dx
    233             d2 = Func(XYZ,Amat)
    234             x += dx
    235             deriv[i][j] = (d1-d2)/(2*dx)
    236     return deriv
     239        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
     240            XYZ[j] += x
     241            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
     242            XYZ[j] -= 2*x
     243            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
     244            XYZ[j] += x
     245            deriv[j][i] = (d1-d2)/(2*dx)
     246    return deriv.flatten()
    237247
    238248def getRestAngle(XYZ,Amat):

  • trunk/GSASIIphsGUI.py

    r817 r820  
    40204020                        UseList[histoName] = {'Histogram':histoName,'Show':False,
    40214021                            'Scale':[1.0,False],'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{}],
    4022                             'Size':['isotropic',[1.,1.,1.0],[False,False,False],[0,0,1],
    4023                                 [4.,4.,4.,0.,0.,0.],6*[False,]],
     4022                            'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1],
     4023                                [1.,1.,1.,0.,0.,0.],6*[False,]],
    40244024                            'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1],
    40254025                                NShkl*[0.01,],NShkl*[False,]],

  • trunk/GSASIIrestrGUI.py

    r819 r820  
    200200                        Lists[listName].append([Ids[x],Types[x],Coords[x],])
    201201        bond = 1.54
    202         dlg = G2phG.SingleFloatDialog(G2frame,'Distance','Enter restraint distance for bond',bond,[1.,4.],'%.4f')
     202        dlg = G2phG.SingleFloatDialog(G2frame,'Distance','Enter restraint distance for bond',bond,[0.01,4.],'%.4f')
    203203        if dlg.ShowModal() == wx.ID_OK:
    204204            bond = dlg.GetValue()
    … …  
    682682           
    683683        def OnUseData(event):
     684            Obj = event.GetEventObject()
    684685            restData['Use'] = Obj.GetValue()
    685686
    … …  
    782783            table = []
    783784            rowLabels = []
     785            bad = []
    784786            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,3',]
    785787            if 'macro' in General['Type']:
    786788                colLabels = ['(res) A - (res) B','calc','obs','esd']
    787789                for i,[indx,ops,obs,esd] in enumerate(bondList):
    788                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
    789                     name = ''
    790                     for atom in atoms:
    791                         name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' - '
    792                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    793                     calc = G2mth.getRestDist(XYZ,Amat)
    794                     table.append([name[:-3],calc,obs,esd])
    795                     rowLabels.append(str(i))               
     790                    try:
     791                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
     792                        name = ''
     793                        for atom in atoms:
     794                            name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' - '
     795                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     796                        calc = G2mth.getRestDist(XYZ,Amat)
     797                        table.append([name[:-3],calc,obs,esd])
     798                        rowLabels.append(str(i))               
     799                    except KeyError:
     800                        print '**** WARNING - missing atom - restraint deleted ****'
     801                        bad.append(i)
    796802            else:
    797803                colLabels = ['A+SymOp - B+SymOp','calc','obs','esd']
    798                 try:
    799                     for i,[indx,ops,obs,esd] in enumerate(bondList):
     804                for i,[indx,ops,obs,esd] in enumerate(bondList):
     805                    try:
    800806                        names = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    801807                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    … …  
    804810                        table.append([names[0]+'+('+ops[0]+') - '+names[1]+'+('+ops[1]+')',calc,obs,esd])
    805811                        rowLabels.append(str(i))
    806                 except ValueError:              #patch - old restraints!
    807                     bondList = []
     812                    except KeyError:
     813                        print '**** WARNING - missing atom - restraint deleted ****'
     814                        bad.append(i)
     815            if len(bad):
     816                bad.reverse()
     817                for ibad in bad:
     818                    del bondList[ibad]
    808819            bondTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    809820            Bonds = G2gd.GSGrid(BondRestr)
    … …  
    880891            table = []
    881892            rowLabels = []
     893            bad = []
    882894            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,2',]
    883895            if 'macro' in General['Type']:
    884896                colLabels = ['(res) A - (res) B - (res) C','calc','obs','esd']
    885897                for i,[indx,ops,obs,esd] in enumerate(angleList):
    886                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
    887                     name = ''
    888                     for atom in atoms:
    889                         name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' - '
    890                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    891                     calc = G2mth.getRestAngle(XYZ,Amat)
    892                     table.append([name[:-3],calc,obs,esd])
    893                     rowLabels.append(str(i))                               
     898                    try:
     899                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
     900                        name = ''
     901                        for atom in atoms:
     902                            name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' - '
     903                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     904                        calc = G2mth.getRestAngle(XYZ,Amat)
     905                        table.append([name[:-3],calc,obs,esd])
     906                        rowLabels.append(str(i))                               
     907                    except KeyError:
     908                        print '**** WARNING - missing atom - restraint deleted ****'
     909                        bad.append(i)
    894910            else:
    895911                colLabels = ['A+SymOp - B+SymOp - C+SymOp','calc','obs','esd']
    896912                for i,[indx,ops,obs,esd] in enumerate(angleList):
    897                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    898                     name = atoms[0]+'+('+ops[0]+') - '+atoms[1]+'+('+ops[1]+') - '+atoms[2]+ \
    899                     '+('+ops[2]+')'
    900                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    901                     XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
    902                     calc = G2mth.getRestAngle(XYZ,Amat)
    903                     table.append([name,calc,obs,esd])
    904                     rowLabels.append(str(i))
     913                    try:
     914                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
     915                        name = atoms[0]+'+('+ops[0]+') - '+atoms[1]+'+('+ops[1]+') - '+atoms[2]+ \
     916                        '+('+ops[2]+')'
     917                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     918                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     919                        calc = G2mth.getRestAngle(XYZ,Amat)
     920                        table.append([name,calc,obs,esd])
     921                        rowLabels.append(str(i))
     922                    except KeyError:
     923                        print '**** WARNING - missing atom - restraint deleted ****'
     924                        bad.append(i)
     925            if len(bad):
     926                bad.reverse()
     927                for ibad in bad:
     928                    del angleList[ibad]
    905929            angleTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    906930            Angles = G2gd.GSGrid(AngleRestr)
    … …  
    968992            table = []
    969993            rowLabels = []
     994            bad = []
    970995            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,2',]
    971996            if 'macro' in General['Type']:
    972997                colLabels = ['(res) atom','calc','obs','esd']
    973998                for i,[indx,ops,obs,esd] in enumerate(planeList):
    974                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
    975                     name = ''
    976                     for a,atom in enumerate(atoms):
    977                         name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' - '
    978                         if (a+1)%3 == 0:
    979                             name += '\n'
    980                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    981                     calc = G2mth.getRestPlane(XYZ,Amat)
    982                     table.append([name[:-3],calc,obs,esd])
    983                     rowLabels.append(str(i))
     999                    try:
     1000                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
     1001                        name = ''
     1002                        for a,atom in enumerate(atoms):
     1003                            name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' - '
     1004                            if (a+1)%3 == 0:
     1005                                name += '\n'
     1006                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     1007                        calc = G2mth.getRestPlane(XYZ,Amat)
     1008                        table.append([name[:-3],calc,obs,esd])
     1009                        rowLabels.append(str(i))
     1010                    except KeyError:
     1011                        print '**** WARNING - missing atom - restraint deleted ****'
     1012                        bad.append(i)
    9841013            else:                               
    9851014                colLabels = ['atom+SymOp','calc','obs','esd']
    9861015                for i,[indx,ops,obs,esd] in enumerate(planeList):
    987                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    988                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    989                     XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
    990                     calc = G2mth.getRestPlane(XYZ,Amat)
    991                     name = ''
    992                     for a,atom in enumerate(atoms):
    993                         name += atom+'+ ('+ops[a]+'),'
    994                         if (a+1)%3 == 0:
    995                             name += '\n'
    996                     table.append([name[:-1],calc,obs,esd])
    997                     rowLabels.append(str(i))
     1016                    try:
     1017                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
     1018                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     1019                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     1020                        calc = G2mth.getRestPlane(XYZ,Amat)
     1021                        name = ''
     1022                        for a,atom in enumerate(atoms):
     1023                            name += atom+'+ ('+ops[a]+'),'
     1024                            if (a+1)%3 == 0:
     1025                                name += '\n'
     1026                        table.append([name[:-1],calc,obs,esd])
     1027                        rowLabels.append(str(i))
     1028                    except KeyError:
     1029                        print '**** WARNING - missing atom - restraint deleted ****'
     1030                        bad.append(i)
     1031            if len(bad):
     1032                bad.reverse()
     1033                for ibad in bad:
     1034                    del planeList[ibad]
    9981035            planeTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    9991036            Planes = G2gd.GSGrid(PlaneRestr)
    … …  
    10711108            table = []
    10721109            rowLabels = []
     1110            bad = []
    10731111            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,2',]
    10741112            if 'macro' in General['Type']:
    10751113                colLabels = ['(res) O (res) A (res) B (res) C','calc','obs','esd']
    10761114                for i,[indx,ops,obs,esd] in enumerate(volumeList):
    1077                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
    1078                     name = ''
    1079                     for atom in atoms:
    1080                         name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' '
    1081                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    1082                     calc = G2mth.getRestChiral(XYZ,Amat)
    1083                     table.append([name,calc,obs,esd])
    1084                     rowLabels.append(str(i))
     1115                    try:
     1116                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
     1117                        name = ''
     1118                        for atom in atoms:
     1119                            name += '('+atom[1]+atom[0].strip()+atom[2]+') '+atom[3]+' '
     1120                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     1121                        calc = G2mth.getRestChiral(XYZ,Amat)
     1122                        table.append([name,calc,obs,esd])
     1123                        rowLabels.append(str(i))
     1124                    except KeyError:
     1125                        print '**** WARNING - missing atom - restraint deleted ****'
     1126                        bad.append(i)
    10851127            else:
    10861128                colLabels = ['O+SymOp  A+SymOp  B+SymOp  C+SymOp)','calc','obs','esd']
    10871129                for i,[indx,ops,obs,esd] in enumerate(volumeList):
    1088                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    1089                     name = atoms[0]+'+('+ops[0]+') '+atoms[1]+'+('+ops[1]+') '+atoms[2]+ \
    1090                         '+('+ops[2]+') '+atoms[3]+'+('+ops[3]+')'
    1091                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    1092                     XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
    1093                     calc = G2mth.getRestChiral(XYZ,Amat)
    1094                     table.append([name,calc,obs,esd])
    1095                     rowLabels.append(str(i))
     1130                    try:
     1131                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
     1132                        name = atoms[0]+'+('+ops[0]+') '+atoms[1]+'+('+ops[1]+') '+atoms[2]+ \
     1133                            '+('+ops[2]+') '+atoms[3]+'+('+ops[3]+')'
     1134                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     1135                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     1136                        calc = G2mth.getRestChiral(XYZ,Amat)
     1137                        table.append([name,calc,obs,esd])
     1138                        rowLabels.append(str(i))
     1139                    except KeyError:
     1140                        print '**** WARNING - missing atom - restraint deleted ****'
     1141                        bad.append(i)
     1142            if len(bad):
     1143                bad.reverse()
     1144                for ibad in bad:
     1145                    del volumeList[ibad]
    10961146            volumeTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    10971147            Volumes = G2gd.GSGrid(ChiralRestr)
    … …  
    11561206            table = []
    11571207            rowLabels = []
     1208            bad = []
    11581209            Types = 2*[wg.GRID_VALUE_STRING,]+4*[wg.GRID_VALUE_FLOAT+':10,2',]
    11591210            if 'macro' in General['Type']:
    11601211                colLabels = ['(res) A  B  C  D','coef name','torsion','obs E','restr','esd']
    11611212                for i,[indx,ops,cofName,esd] in enumerate(torsionList):
    1162                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
    1163                     name = '('+atoms[2][1]+atoms[2][0].strip()+atoms[2][2]+')'
    1164                     for atom in atoms:
    1165                         name += '  '+atom[3]
    1166                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    1167                     tor = G2mth.getRestTorsion(XYZ,Amat)
    1168                     restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
    1169                     table.append([name,cofName,tor,calc,restr,esd])
    1170                     rowLabels.append(str(i))
     1213                    try:
     1214                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
     1215                        name = '('+atoms[2][1]+atoms[2][0].strip()+atoms[2][2]+')'
     1216                        for atom in atoms:
     1217                            name += '  '+atom[3]
     1218                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     1219                        tor = G2mth.getRestTorsion(XYZ,Amat)
     1220                        restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
     1221                        table.append([name,cofName,tor,calc,restr,esd])
     1222                        rowLabels.append(str(i))
     1223                    except KeyError:
     1224                        print '**** WARNING - missing atom - restraint deleted ****'
     1225                        bad.append(i)
     1226            if len(bad):
     1227                bad.reverse()
     1228                for ibad in bad:
     1229                    del torsionList[ibad]
    11711230            torsionTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    11721231            Torsions = G2gd.GSGrid(TorsionRestr)
    … …  
    12491308            table = []
    12501309            rowLabels = []
     1310            bad = []
    12511311            Types = 2*[wg.GRID_VALUE_STRING,]+5*[wg.GRID_VALUE_FLOAT+':10,2',]
    12521312            if 'macro' in General['Type']:
    12531313                colLabels = ['(res) A  B  C  D  E','coef name','phi','psi','obs E','restr','esd']
    12541314                for i,[indx,ops,cofName,esd] in enumerate(ramaList):
    1255                     atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
    1256                     name = '('+atoms[3][1]+atoms[3][0].strip()+atoms[3][2]+')'
    1257                     for atom in atoms:
    1258                         name += '  '+atom[3]
    1259                     XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
    1260                     phi,psi = G2mth.getRestRama(XYZ,Amat)
    1261                     restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
    1262                     table.append([name,cofName,phi,psi,calc,restr,esd])
    1263                     rowLabels.append(str(i))
     1315                    try:
     1316                        atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,0,4)
     1317                        name = '('+atoms[3][1]+atoms[3][0].strip()+atoms[3][2]+')'
     1318                        for atom in atoms:
     1319                            name += '  '+atom[3]
     1320                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cx,3))
     1321                        phi,psi = G2mth.getRestRama(XYZ,Amat)
     1322                        restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
     1323                        table.append([name,cofName,phi,psi,calc,restr,esd])
     1324                        rowLabels.append(str(i))
     1325                    except KeyError:
     1326                        print '**** WARNING - missing atom - restraint deleted ****'
     1327                        bad.append(i)
     1328            if len(bad):
     1329                bad.reverse()
     1330                for ibad in bad:
     1331                    del ramaList[ibad]
    12641332            ramaTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    12651333            Ramas = G2gd.GSGrid(RamaRestr)

  • trunk/GSASIIstruct.py

    r819 r820  
    293293                if hist not in Histograms and Phase['Histograms'][hist]['Use']:
    294294                    Histograms[hist] = allHistograms[hist]
    295                     #future restraint, etc. histograms here           
    296295                    hId = histoList.index(hist)
    297296                    Histograms[hist]['hId'] = hId
    … …  
    512511################################################################################       
    513512                   
    514 def GetPhaseData(PhaseData,RestraintDict=None,Print=True,pFile=None):
     513def GetPhaseData(PhaseData,RestraintDict={},Print=True,pFile=None):
    515514           
    516515    def PrintFFtable(FFtable):
    … …  
    600599        print >>pFile,ptstr
    601600       
    602     def PrintRestraints(phaseRest):
    603         if phaseRest:
    604             print >>pFile,'\n Restraints:'
    605             names = ['Bonds','Angles','Planes','Volumes','Torsions','Ramas']
    606             for i,name in enumerate(['Bond','Angle','Plane','Chiral','Torsion','Rama']):
    607                 itemRest = phaseRest[name]
    608                 if itemRest[names[i]]:
    609                     print >>pFile,'\n  %30s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
    610                     print >>pFile,'  atoms(symOp), calc, obs, sig: '
    611                     for item in phaseRest[name][names[i]]:
    612                         print item
    613 #                        text = '   '
    614 #                        for a,at in enumerate(item[0]):
    615 #                            text += at+'+('+item[1][a]+') '
    616 #                            if (a+1)%5 == 0:
    617 #                                text += '\n   '
    618 #                        print >>pFile,text,' %.3f %.3f %.3f'%(item[3],item[4],item[5])
    619        
    620601    if Print:print  >>pFile,' Phases:'
    621602    phaseVary = []
    … …  
    640621        BLtables.update(BLtable)
    641622        Atoms = PhaseData[name]['Atoms']
     623        AtLookup = G2mth.FillAtomLookUp(Atoms)
    642624        try:
    643625            PawleyRef = PhaseData[name]['Pawley ref']
    … …  
    648630        cell = General['Cell']
    649631        A = G2lat.cell2A(cell[1:7])
     632        Amat,Bmat = G2lat.cell2AB(cell[1:7])
    650633        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
    651634        if cell[0]:
    … …  
    727710                PrintTexture(textureData)
    728711                if name in RestraintDict:
    729                     PrintRestraints(RestraintDict[name])
     712                    PrintRestraints(Amat,SGData,General['AtomPtrs'],Atoms,AtLookup,RestraintDict[name],pFile)
    730713                   
    731714        elif PawleyRef:
    … …  
    781764        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
    782765    return A,sigA
     766       
     767def PrintRestraints(Amat,SGData,AtPtrs,Atoms,AtLookup,phaseRest,pFile):
     768    if phaseRest:
     769        cx,ct = AtPtrs[:2]
     770        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
     771            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas']]
     772        for name,rest in names:
     773            itemRest = phaseRest[name]
     774            if itemRest[rest] and itemRest['Use']:
     775                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
     776                if name in ['Bond','Angle','Plane','Chiral']:
     777                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
     778                    for indx,ops,obs,esd in itemRest[rest]:
     779                        try:
     780                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
     781                            AtName = ''
     782                            for i,Aname in enumerate(AtNames):
     783                                AtName += Aname
     784                                if ops[i] == '1':
     785                                    AtName += '-'
     786                                else:
     787                                    AtName += '+('+ops[i]+')-'
     788                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
     789                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     790                            if name == 'Bond':
     791                                calc = G2mth.getRestDist(XYZ,Amat)
     792                            elif name == 'Angle':
     793                                calc = G2mth.getRestAngle(XYZ,Amat)
     794                            elif name == 'Plane':
     795                                calc = G2mth.getRestPlane(XYZ,Amat)
     796                            elif name == 'Chiral':
     797                                calc = G2mth.getRestChiral(XYZ,Amat)
     798                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
     799                        except KeyError:
     800                            del itemRest[rest]
     801                else:
     802                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
     803                    coeffDict = itemRest['Coeff']
     804                    for indx,ops,cofName,esd in enumerate(torsionList):
     805                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
     806                        AtName = ''
     807                        for i,Aname in enumerate(AtNames):
     808                            AtName += Aname+'+('+ops[i]+')-'
     809                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
     810                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     811                        if name == 'Torsion':
     812                            tor = G2mth.getRestTorsion(XYZ,Amat)
     813                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
     814                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(AtName[:-1],calc,obs,esd,(obs-calc)/esd,tor)
     815                        else:
     816                            phi,psi = G2mth.getRestRama(XYZ,Amat)
     817                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
     818                            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)
    783819       
    784820def getCellEsd(pfx,SGData,A,covData):
    … …  
    840876    return cellSig           
    841877   
    842 def SetPhaseData(parmDict,sigDict,Phases,covData,pFile=None):
     878def SetPhaseData(parmDict,sigDict,Phases,covData,RestraintDict=None,pFile=None):
    843879   
    844880    def PrintAtomsAndSig(General,Atoms,atomsSig):
    … …  
    929965        SGData = General['SGData']
    930966        Atoms = Phase['Atoms']
     967        AtLookup = G2mth.FillAtomLookUp(Atoms)
    931968        cell = General['Cell']
     969        Amat,Bmat = G2lat.cell2AB(cell[1:7])
    932970        pId = Phase['pId']
    933971        pfx = str(pId)+'::'
    … …  
    9821020        else:
    9831021            atomsSig = {}
    984             if General['Type'] == 'nuclear':
     1022            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
    9851023                for i,at in enumerate(Atoms):
    9861024                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
    … …  
    10221060                    SHtextureSig[name] = sigDict[aname]
    10231061            PrintSHtextureAndSig(textureData,SHtextureSig)
     1062        if phase in RestraintDict:
     1063            PrintRestraints(Amat,SGData,General['AtomPtrs'],Atoms,AtLookup,RestraintDict[phase],pFile)
     1064                   
    10241065
    10251066################################################################################
    … …  
    19581999################################################################################
    19592000
    1960 def penaltyFxn(Phases,parmDict,varyList):
     2001def penaltyFxn(HistoPhases,parmDict,varyList):
     2002    Histograms,Phases,restraintDict = HistoPhases
    19612003    pNames = []
    19622004    pVals = []
     2005    pWt = []
    19632006    negWt = {}
    19642007    for phase in Phases:
    1965         negWt[Phases[phase]['pId']] = Phases[phase]['General']['Pawley neg wt']
     2008        pId = Phases[phase]['pId']
     2009        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
     2010        General = Phases[phase]['General']
     2011        SGData = General['SGData']
     2012        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
     2013        Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
     2014        if phase not in restraintDict:
     2015            continue
     2016        phaseRest = restraintDict[phase]
     2017        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
     2018            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas']]
     2019        for name,rest in names:
     2020            itemRest = phaseRest[name]
     2021            if itemRest[rest] and itemRest['Use']:
     2022                wt = itemRest['wtFactor']
     2023                if name in ['Bond','Angle','Plane','Chiral']:
     2024                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
     2025                        pNames.append(str(pId)+':'+name+':'+str(i))
     2026                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
     2027                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     2028                        if name == 'Bond':
     2029                            calc = G2mth.getRestDist(XYZ,Amat)
     2030                        elif name == 'Angle':
     2031                            calc = G2mth.getRestAngle(XYZ,Amat)
     2032                        elif name == 'Plane':
     2033                            calc = G2mth.getRestPlane(XYZ,Amat)
     2034                        elif name == 'Chiral':
     2035                            calc = G2mth.getRestChiral(XYZ,Amat)
     2036                        pVals.append(obs-calc)
     2037                        pWt.append(wt/esd**2)
     2038                else:
     2039                    coeffDict = itemRest['Coeff']
     2040                    for i,[indx,ops,cofName,esd] in enumerate(torsionList):
     2041                        pNames.append(str(pId)+':'+name+':'+str(i))
     2042                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
     2043                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
     2044                        if name == 'Torsion':
     2045                            tor = G2mth.getRestTorsion(XYZ,Amat)
     2046                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
     2047                        else:
     2048                            phi,psi = G2mth.getRestRama(XYZ,Amat)
     2049                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
     2050                        pVals.append(obs-calc)
     2051                        pWt.append(wt/esd**2)
     2052         
    19662053    for item in varyList:
    19672054        if 'PWLref' in item and parmDict[item] < 0.:
    19682055            pId = int(item.split(':')[0])
    19692056            pNames.append(item)
    1970 #            pVals.append(negWt[pId]*parmDict[item]**2)
    1971             pVals.append(-negWt[pId]*parmDict[item])
     2057            pVals.append(-parmDict[item])
     2058            pWt.append(negWt[pId])
    19722059    pVals = np.array(pVals)
    1973     return pNames,pVals
    1974    
    1975 def penaltyDeriv(pNames,pVal,Phases,parmDict,varyList):
     2060    pWt = np.array(pWt)
     2061    return pNames,pVals,pWt
     2062   
     2063def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
     2064    Histograms,Phases,restraintDict = HistoPhases
    19762065    pDerv = np.zeros((len(varyList),len(pVal)))
    1977     negWt = {}
    19782066    for phase in Phases:
    1979         negWt[Phases[phase]['pId']] = np.sqrt(Phases[phase]['General']['Pawley neg wt'])
     2067        pId = Phases[phase]['pId']
     2068        General = Phases[phase]['General']
     2069        SGData = General['SGData']
     2070        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
     2071        Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
     2072        phaseRest = restraintDict[phase]
     2073        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
     2074            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas'}
     2075        for ip,pName in enumerate(pNames):
     2076            pnames = pName.split(':')
     2077            if pId == int(pnames[0]):
     2078                name = pnames[1]
     2079                id = int(pnames[2])
     2080                itemRest = phaseRest[name]
     2081                if name in ['Bond','Angle','Plane','Chiral']:
     2082                    indx,ops,obs,esd = itemRest[names[name]][id]
     2083                    dNames = []
     2084                    for ind in indx:
     2085                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
     2086                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
     2087                    if name == 'Bond':
     2088                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
     2089                    elif name == 'Angle':
     2090                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
     2091                    elif name == 'Plane':
     2092                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
     2093                    elif name == 'Chiral':
     2094                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
     2095                else:
     2096                    coeffDict = itemRest['Coeff']
     2097                    for indx,ops,obs,esd in itemRest[id]:
     2098                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
     2099#                        if name == 'Torsion':
     2100#                            tor = G2mth.getRestTorsion(XYZ,Amat)
     2101#                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
     2102#                        else:
     2103#                            phi,psi = G2mth.getRestRama(XYZ,Amat)
     2104#                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
     2105                for dName,drv in zip(dNames,deriv):
     2106                    try:
     2107                        ind = varyList.index(dName)
     2108                        pDerv[ind][ip] += drv
     2109                    except ValueError:
     2110                        pass
    19802111    for i,item in enumerate(varyList):
    19812112        if item in pNames:
    19822113            pId = int(item.split(':')[0])
    1983 #            pDerv[i][pNames.index(item)] += 2.*negWt[pId]*parmDict[item]
    1984             pDerv[i][pNames.index(item)] += negWt[pId]
     2114            pDerv[i][pNames.index(item)] += 1.
    19852115    return pDerv
    19862116
    … …  
    30723202            dMdv = dMdvh
    30733203           
    3074     pNames,pVals = penaltyFxn(Phases,parmdict,varylist)
    3075     if np.sum(pVals):
    3076         dpdv = penaltyDeriv(pNames,pVals,Phases,parmdict,varylist)
     3204    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
     3205    if np.any(pVals):
     3206        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
    30773207        dMdv = np.concatenate((dMdv.T,dpdv.T)).T
     3208       
    30783209    return dMdv
    30793210
    … …  
    31613292        else:
    31623293            continue        #skip non-histogram entries
    3163     pNames,pVals = penaltyFxn(Phases,parmdict,varylist)
    3164     if np.sum(pVals):
    3165         dpdv = penaltyDeriv(pNames,pVals,Phases,parmdict,varylist)
    3166         Vec += np.sum(dpdv,axis=1)
    3167         Hess += np.inner(dpdv,dpdv)
     3294    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
     3295    if np.any(pVals):
     3296        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
     3297        Vec += np.sum(dpdv*pWt*pVals,axis=1)
     3298        Hess += np.inner(dpdv*pWt,dpdv)
    31683299    return Vec,Hess
    31693300
    … …  
    32763407            dlg.Destroy()
    32773408            raise Exception         #Abort!!
    3278     pDict,pVals = penaltyFxn(Phases,parmdict,varylist)
    3279     if np.sum(pVals):
    3280         print 'Penalty function :',np.sum(pVals),' on ',len(pVals),' terms'
    3281         M = np.concatenate((M,pVals))
     3409    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
     3410    if np.any(pVals):
     3411        print 'Penalty function: %.3f on %d terms'%(np.sum(pWt*pVals**2),len(pVals))
     3412        M = np.concatenate((M,pWt*pVals))
    32823413    return M
    32833414                       
    … …  
    34173548    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
    34183549    G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
    3419     SetPhaseData(parmDict,sigDict,Phases,covData,printFile)
     3550    SetPhaseData(parmDict,sigDict,Phases,covData,restraintDict,printFile)
    34203551    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
    34213552    SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
Note: See TracChangeset for help on using the changeset viewer.