[gmx-users] Re: Re: Positive potential energy for TFE solvent

Mark Abraham Mark.Abraham at anu.edu.au
Thu Nov 17 06:48:35 CET 2011


On 16/11/2011 5:02 PM, Harpreet Basra wrote:
> Hi Mark,
> Sorry my last mail was incomplete...here is the complete one!
>
>     On 16/11/2011 1:18 AM, Harpreet Basra wrote:
>     > Hi Mark,
>     >
>     > Thanks for the quick reply. But i have already done what u
>     suggested.
>     >
>     >
>     >
>     >     On 15/11/2011 6:06 PM, Harpreet Basra wrote:
>     > > Hi
>     > > I am still stuck with same problem of obtaining positive potential
>     > > energy.
>     > > >>On 11/11/2011 5:07 PM, Harpreet Basra wrote:
>     > > >> Hi
>     > > >>
>     > > >> I am trying to generate an equilibrated box of 216 TFE
>     >     molecules.To
>     > > >> generate the 216 TFE molecule box i performed following steps:
>     > > >
>     > > >A suggested workflow can be found here
>     > > >http://www.gromacs.org/Documentation/How-tos/Non-Water_Solvation
>     > > I have been following this link only.
>     > > >
>     > > >
>     > > >> 1) I got the tfe.gro file and created a cubic box of edge
>     >     length =
>     > > >> 0.516 nm containing 1 TFE molecule (at its center), using the
>     > > >> following command:
>     > > >>
>     > > >>>> editconf -f tfe.gro -c -o tfe_box.gro -bt cubic -box 0.516
>     > > >>  I chose this length because in the tfe.gro file dimensions
>     >     of the TFE
>     > > >> molecule are 0.516 0.516 0.516.
>     > > >
>     > > >That's not a good reason. Choose a volume and shape that makes
>     >     sense for
>     > > >your target density. Cubic probably doesn't make sense when a
>     > > >rectangular shape is possible. Then you'll probably want to
>     >     choose -nbox
>     > > >differently later.
>     > > I chose a rectangular box too. still i get a positive value for
>     >     PE and
>     > > moreover all the molecules move towards two opposite walls of
>     >     the box.
>     > > I am not sure that the way I am using the genconf command is the
>     > > correct way. because I have tried every other possibility for not
>     > > getting a positive potential, with no success. So here are my .gro
>     > > file and the topology file for TFE.
>     > > *****tfe.gro file*****
>     > > 7
>     > >
>     > > 1TFE   F1T   1   0.444   0.344   0.246
>     > >
>     > > 1TFE   CT     2   0.334  0.245   0.246
>     > >
>     > > 1TFE   F2T    3   0.350  0.160   0.364
>     > >
>     > > 1TFE   F3T    4   0.350  0.160   0.127
>     > >
>     > > 1TFE   CH2T  5  0.187  0.326   0.246
>     > >
>     > > 1TFE   OT      6  0.075  0.220   0.246
>     > >
>     > > 1TFE   HT      7  -0.019 0.266   0.246
>     > >
>     > > 0.49174   0.49174   0.49174
>     > >
>     > > ****topology file****
>     > >
>     > > [ moleculetype ]
>     > >
>     > > ; Name nrexcl
>     > >
>     > > TFE 3
>     > >
>     > > [ atoms ]
>     > >
>     > > ; nr type resnr resid atom cgnr charge mass
>     > >
>     > > 1 FTFE    1  TFE  F1T   1   -0.170   18.9984
>     > >
>     > > 2 CTFE    1  TFE  CT    1    0.452   12.0110
>     > >
>     > > 3 FTFE    1  TFE  F2T   1   -0.170   18.9984
>     > >
>     > > 4 FTFE    1  TFE  F3T   1   -0.170   18.9984
>     > >
>     > > 5 CHTFE 1  TFE  CH2T 1    0.273   14.0270
>     > >
>     > > 6 OTFE   1  TFE  OT     1   -0.625   15.9994
>     > >
>     > > 7 H          1  TFE  HT     1   0.410   1.0080
>     > >
>     > > [ bonds ]
>     > >
>     > > ; ai aj fu c0, c1, ...
>     > >
>     > > 2 1 2 0.133 3380866.9 0.133 3380866.9 ; C1 F1
>     > >
>     > > 2 3 2 0.133 3380866.9 0.133 3380866.9 ; C1 F2
>     > >
>     > > 2 4 2 0.133 3380866.9 0.133 3380866.9 ; C1 F3
>     > >
>     > > 2 5 2 0.153 7150000.0 0.153 7150000.0 ; C1 C2
>     > >
>     > > 5 6 2 0.143 8180000.0 0.143 8180000.0 ; C2 O
>     > >
>     > > 6 7 2 0.100 15700000.0 0.100 15700000.0 ; O H
>     > >
>     > > [ pairs ]
>     > >
>     > > ; ai aj fu c0, c1, ...
>     > >
>     > > 1 6 1 ; F1 O
>     > >
>     > > 2 7 1 ; C1 H
>     > >
>     > > 3 6 1 ; F2 O
>     > >
>     > > 4 6 1 ; F3 O
>     > >
>     > > [ angles ]
>     > >
>     > > ; ai aj ak fu c0, c1, ...
>     > >
>     > > 1 2 3 2 109.5 520.0 109.5 520.0 ; F1 C1 F2
>     > >
>     > > 1 2 4 2 109.5 520.0 109.5 520.0 ; F1 C1 F3
>     > >
>     > > 1 2 5 2 109.5 520.0 109.5 520.0 ; F1 C1 C2
>     > >
>     > > 3 2 4 2 109.5 520.0 109.5 520.0 ; F2 C1 F3
>     > >
>     > > 3 2 5 2 109.5 520.0 109.5 520.0 ; F2 C1 C2
>     > >
>     > > 4 2 5 2 109.5 520.0 109.5 520.0 ; F3 C1 C2
>     > >
>     > > 2 5 6 2 109.5 520.0 109.5 520.0 ; C1 C2 O
>     > >
>     > > 5 6 7 2 109.5 450.0 109.5 450.0 ; C2 O H
>     > >
>     > > [ dihedrals ]
>     > >
>     > > ; ai aj ak al fu c0, c1, m, ...
>     > >
>     > > 6 5 2 1 1 0.0 5.9 3 0.0 5.9 3 ; dih O C2 C1 F1
>     > >
>     > > 2 5 6 7 1 0.0 1.3 3 0.0 1.3 3 ; dih C1 C2 O H
>     > >
>     > > and to construct a box of TFE solvent i took the tfe.gro file and
>     > > replicated the TFE molecule by using
>     > > genconf -f tfe.gro -o tfe_sol.gro -rot -nbox 6
>     > > can u plz suggest is it that I am using genconf in a wrong way
>     >     that it
>     > > is causing this problem? I am not sure how many molecules (-nbox
>     > > option in genconf) should i keep in the box in order to get a mass
>     > > density of 1383g/L for TFE.
>     >
>     >     That link says "Work out how much volume a single molecule would
>     >     have in
>     >     the box of your chosen density and size. Useeditconf
>     > <http://www.gromacs.org/editconf>to place a box of that size
>     >     around your
>     >     single molecule." It does not seem to me that you have done
>     this.
>     >
>     >     Mark
>     >
>     >
>     > I did place the *single molecule* in a box of size required to get a
>     > density of 1383 g/L. I also checked the density of the solvent box
>     > (containing 216 molecules after NVT equilibration for 200 ps) I
>     > constructed the average value comes out to be 1397 g/L with a std
>     > deviation of 30 g/L, thus it seems range. Moreover, the potential
>     > energy of this one molecule (tfe.gro) was coming out to be
>     > highly negative (-6.4E+08 kJ/mol). But on generating a solvent
>     system
>     > with 216 TFE molecules the energy becomes (1.9E+04 kJ/mol).
>     >
>
>     Since a single molecule looks OK, the usual reasons for a large
>     positive
>     PE would be atomic overlap (visualise the result of genconf -
>     perhaps a
>     cubic box is not a good shape), or invalid inter-molecular non-bonded
>     parameters (place two boxes 0.5nm apart and see what the energy does).
>
>     Mark
>
> I tried genconf taking (1) a rectangular box and (2) keeping a 
> distance of 0.5nm between two boxes. What i observe in both the cases 
> is that all the molecules gather at two opposite walls of the box 
> creating a large empty space at the center. The reason for this seems 
> to be formation of hydrogen bonds between the two separate clusters of 
> TFE on both walls. Since on calculating the hydrogen bonds between the 
> TFEmolecules the average number comes out to be 189 +/- 34. Also on 
> visualizing the system in VMD i see a strong network of hydrogen bonds 
> in the two separated clusters.
>
>
> Apart from this u asked me to see the result after genconf. So i 
> calculated the energies just after genconf for the system having a 
> distance of 0.5 nm between two adjacent boxes. again the potential is 
> positive. here are the results afte genconf and energy minimization:
> Energies (kJ/mol)
>         G96Bond       G96Angle        Proper Dih.          
> LJ-14         Coulomb-14
>     2.12101e+02    6.94824e+02    3.88264e-01    6.38056e+01    
> 4.67073e+04
>         LJ (SR)       Coulomb (SR)    Coul. recip.      
> Potential        Pressure (bar)
>    -4.10111e+01   -1.68897e+03   -1.39086e+04    3.20399e+04   
> -1.31929e+02
> so just after genconf there is no atomic overlap since the molecules 
> are far apart,
> the potential is still positive.

Coulomb-14 is far too large. Look at the energy of the pre-miminization 
structure too. If you visualize the result of genconf with the periodic 
box shown (like I've suggested at least once), then you can start to 
address why this intra-molecular contribution is out of control.

Mark

> Thinking of this i tried to visualize the system previously generated 
> by me (the one i reported to u earlier). here the system looks 
> uniformally distributed with a strong network of hydrogen bonds. I 
> confirmed this by calculatein the no. of hydrogen bonds in my 
> equilibrated (NVT) trajectory which comes out to be 199 +/- 11. 
> moreover the different energy terms from the log file are as follows:
> <====  ###############  ==>
> <====  A V E R A G E S  ====>
> <==  ###############  ======>
> Energies (kJ/mol)
>        G96Angle    Proper Dih.          LJ-14     Coulomb-14        LJ 
> (SR)
>     2.73173e+03    4.51905e+02    1.90037e+02    4.85110e+04   
> -3.35996e+03
>    Coulomb (SR)   Coul. recip.      Potential    Kinetic En.   Total 
> Energy
>    -9.37385e+03   -1.52068e+04    2.39441e+04    4.03654e+03    
> 2.79807e+04
>     Temperature Pressure (bar)  Cons. rmsd ()
>     2.99958e+02    2.52412e+02    0.00000e+00
> so do u think that getting a positive total potential energy is a 
> result of some atomic overlap or there are actually any repulsive 
> forces operative in my TFE solvent system? since even after simple 
> genconf when no overlap is possible i get a positive energy.
> Thanks
> Harpreet
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111117/a75bd20c/attachment.html>


More information about the gromacs.org_gmx-users mailing list