Hi Aimee,
1. If I see correctly, in your script you getting the scale factor by
dividing the frame's supercell dimensions by the experimental unit cell
dimensions. This is not good. As a check, print out the scale factors
you get. If your simulation volume is not deviating by more than 0.3%
from simulation, as should be the case, then your scale factors should
all be close to 1. Also visualize the structures after scaling: they
should look "normal".
2. I don't have the script used for fitting the monomers in the 2010
paper (I didn't work on that). Offline I'll send you the email of the
first author of that paper.
3. Regarding the method for calculating B-factors, well it depends on
what you want to get. The "reverse symmetry" B-factors in Janowski et
al., J. Am. Chem. Soc. 2013, 135, 7938-7948 were calculated by method 1.
I assume that that's also how the "lattice b-factors" were done in the
2010 paper.
Pawel
On 11/08/2013 04:19 AM, Yongxiu Li wrote:
> hello,Pawel,
> I am sorry to trouble you again!Thank you for your attention!
> In my work, the volume rescaling was anisotropic so I think I need
> three scale factor of x,y, z vector . I change the script to
> RevSym_mod.py (see attachment ) .After that , the rms fluctuation of
> a-carbons for asu and lattice decrease so much at once (see figure 1)
> but unlucky they are so different from the result of 2010 Paper and
> the asu_rmsd and lattice_rmsd become about 7.0 ? which is so high. So
> I want to know how to do the least squared minimization problem which
> 2010 Paper used as you said. Can you give me this script or tell me
> how to gain it?
> Then I am confused to compute the bfactor. One method is trajin all
> 108 monomers trajectory and compute the bfactor (see method1) and the
> second method is that trajin in each monomer trajectory and compute
> its bfactor then average these 108 bfactors . We can see that the
> lattice_bfactor from method2(see figure3 ) is smaller than that (see
> figure 2) from method 1 but asu_bfator is bigger . So I want to know
> which one is right ? (BTW, before gaining the figure 2 and figure 3,
> we didn't do the rescaling )
> Thank you very much!
> Best wishes!
> Aimee Li
>
>
> method1 :
> #################################
> # Calculate RevSymm B-factors #
> #################################
> f=open('ctraj_bfactor_lat','w')
> f.write('parm %s\n' %topo)
> for i in range(unitcells):
> for j in range(asymunits):
> f.write('trajin RevSym_%02d_%02d.nc\n' %(i+1,j+1))
> f.write('atomicfluct out bfac_lat_calpha.dat :1-64.CA byatom bfactor\n' )
> f.close()
> os.system('cpptraj <ctraj_bfactor_lat')
>
> method2:
> #################################
> # Calculate RotFit B-factors #
> #################################
> for i in range(unitcells):
> for j in range(asymunits):
> f=open('ctraj_bfactor_asu-bk_%02d_%02d' %(i+1,j+1),'w')
> f.write('parm %s\n' %topo)
> f.write('trajin ../RevSym_%02d_%02d.nc\n' %(i+1,j+1))
> f.write('reference AvgCoord_asu.rst7.1\n')
> f.write('rms reference mass :1-64.CA,C,N\n')
> f.write('atomicfluct out bfac_asu_calpha-bk_%02d_%02d.dat :1-64.CA
> byatom bfactor\n' %(i+1,j+1))
> f.close()
> os.system('cpptraj <ctraj_bfactor_asu-bk_%02d_%02d' %(i+1,j+1))
>
> ===============
> On 10/30/2013 09:12 PM, Pawel wrote:
>> Hi Aimee,
>>
>> To get the scaling properly you need something like this:
>>
>> A=40.000 #set this to the experimental value of your
>> supercell's a-vector
>> for frame in range(frames):
>> a=celllen[frame,0] #get the a-vector of the current frame
>> scale=a/A #scale factor (pressure scaling is isotropic in
>> Amber
>> so your scale factor will be the same in all dimensions
>> currMoveVector=MoveVector*scale # scale your MoveVector
>> coords[frame,:,:]=coords[frame,:,:]-MoveVector
>>
>> Try that and see what happens although if you are monitoring your volume
>> and it is not deviating by more than 0.3% or so from the experimental,
>> this should really not affect your B-factors all that much.
>>
>> Besides that I'm really not sure what the reason could be. Possibly the
>> algorithm for fitting the unit cells (the least squares problem for
>> finding the correct origin that I mentioned in previous e-mail) could
>> have made a difference. Possibly your system is sampling an alternate
>> ensemble that is producing different B-factors.
>>
>> Pawel
>>
>>
>>
>> On 10/29/2013 04:03 AM, Yongxiu Li wrote:
>>> hello,Pawel
>>> Thank you very much! I am so sorry for that I read the
>>> RevSym_netcdf.py again and I think I have done the rescaling. I think
>>> these lines do the all rescaling "
>>> for frame in range(frames):
>>> #get frame box
>>> box=celllen[frame,:]
>>> #convert to unit cell box
>>> box=box/[ix,iy,iz]
>>> #add angle information as provided by user
>>> UCbox=hstack((box,SCBox[3:]))
>>> #calc orthogonalization matrix
>>> u,invu=CompXfrm(UCbox)
>>> #matrix product
>>> MoveVector=dot(invu,FracVector).astype(float32)
>>> coords[frame,:,:]=coords[frame,:,:]-MoveVector
>>> # reverse symmetry operate to original asym unit by
>>> # applying the symmetry translation and rotation
>>> coords[:,:,:]=dot( (coords[:,:,:]-t),linalg.inv(s) ) "
>>> I'm not proficient in Python, so whether I understand is right? If
>>> not, can you tell me how to do in python script? If I am right, I
>>> still can't understand that compared with figure8 in Case's paper,why
>>> the bfac_lat_calpha-arf.dat which I gained is so big.
>>> btw,in the attachment is the lattice atomic root mean squared
>>> fluctuations of a-carbons by amber ff99SB-tip3p.
>>> Thank you very much!
>>> Aimee Li
>>>
>>> Atomicrootmeansquared(rms)fluctuationsofa-carbons.
>>>
>>> On 10/16/2013 07:26 AM, Pawel wrote:
>>>> Hi Aimee,
>>>> On 10/15/2013 09:09 AM, Yongxiu Li wrote:
>>>>> 1. if the "bfac_lat_calpha" in figure8 are calculated using a
>>>>> different
>>>>> approach than the one used in BasicAnalysis, how can I do it?
>>>> If I understand correctly (I did not work on that paper myself), the
>>>> alignment of asymmetric units was done by solving a least squares
>>>> minimization problem to find (for each frame) the best location of the
>>>> crystallographic origin in space so that after performing the
>>>> appropriate symmetry operations and translations in reverse on each of
>>>> the asymmetric units, the resulting RMSD between the asymmetric units
>>>> would be minimized. Unfortunately I don't have the code for that
>>>> approach. If it would be very useful to you, let me know and I'll
>>>> investigate.
>>>>
>>>>> So as you said, I need to rescal coordinates.
>>>>> However, after reading this paper I am confused that wether only
>>>>> computing the lattice property (for example:rmsd_latt and
>>>>> bfactor_latt)need rescal coordinates but computing the other
>>>>> property(like rmsd_ASU ,bfactor_ASU and Distance deviation matrices)
>>>>> needn't rescale the coordinates?
>>>> If you are simulating an NPT ensemble, you should rescale the
>>>> coordinates for all of these analyses. This should be done in
>>>> RevSym_netcdf.py though I think the version of that script that you
>>>> have
>>>> does not do so. Like I said these scripts Analysis directory are all
>>>> beta versions. However you should be able to do the rescaling
>>>> easily by
>>>> adding one or two line around line 95 of the script. If you have
>>>> problems, let me know and I'll dig up a version that does this.
>>>>
>>>> Pawel
>>>>
>>>> _______________________________________________
>>>> AMBER mailing list
>>>> AMBER.ambermd.org
>>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>>
>>>
>>>
>>> _______________________________________________
>>> AMBER mailing list
>>> AMBER.ambermd.org
>>> http://lists.ambermd.org/mailman/listinfo/amber
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Nov 08 2013 - 06:30:03 PST
