Changeset 820


Ignore:
Timestamp:
Jan 3, 2013 4:03:22 PM (9 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.