[gmx-users] speed of sound in liquid
Dear all, Does anyone know how to obtain speed of sound, c, in a liquid? c is defined as square root of dP/drho (pressure change with density) at constant entropy My idea is to run a series of NPT simulations to obtain density at each P, to get dP/drho. However, I am not sure how to ensure a constant entropy in the system... Could anyone comment on this? Thanks a lot in advance, E. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] surface tension-crash at step 0
Hi all, I have problem ruunning NVT with Z direction extended to get surface tension of a polymer packed in a cell. I tried the same procedure for an alkane and NVT works well however for the polymer the run crashes at the very first step.. Doe anyone have clue what wrong could be? Many thanks. starting mdrun 'Polymer' 1000 steps, 1.0 ps. step 0 [node6:14338] *** Process received signal *** [node6:14339] *** Process received signal *** [node6:14338] Signal: Segmentation fault (11) [node6:14338] Signal code: Address not mapped (1) [node6:14338] Failing at address: 0x2fbbaca0 [node6:14339] Signal: Segmentation fault (11) [node6:14339] Signal code: Address not mapped (1) [node6:14339] Failing at address: 0x2fbbd2e0 [node6:14339] [ 0] /lib64/libpthread.so.0 [0x3112e0ebe0] [node6:14339] [ 1] /usr/local/gromacs/lib/libgmx_mpi.so.6 [0x2b2fd8203c54] [node6:14339] *** End of error message *** [node6:14338] [ 0] /lib64/libpthread.so.0 [0x3112e0ebe0] [node6:14338] [ 1] /usr/local/gromacs/lib/libgmx_mpi.so.6 [0x2b4b10e86c65] [node6:14338] *** End of error message *** -- mpirun noticed that process rank 3 with PID 14338 on node node6 exited on signal 11 (Segmentation fault). On 4 April 2013 12:59, André Farias de Moura mo...@ufscar.br wrote: just that simple: as soon as you expand the box in the z direction the systems releases the excess pressure by the expansion of the liquid into the evacuated region. On Thu, Apr 4, 2013 at 1:20 PM, Elisabeth katesed...@gmail.com wrote: Dear Dr. Moura, Thank you for your answer. I equilibrated the cell under NPT and then extended the Z direction to get a surface. Since I run NVT and Z is extended, I see no pressure dependence and surface tension values are similar for all P from 50 to 1000 bar. which is attributed to the fact that Z is extended and system tends to expand regardless of the pressure that was imposed in the initial NPT runs. I have no gas, it is liqiuid-vaccum and wanted to know if there is anyway to capture effect of pressure. From your answer I think there is no such possibility... Please comment if the answer is positive, Thank you On 4 April 2013 10:32, André Farias de Moura mo...@ufscar.br wrote: There's no simple answer for that. If you apply a lateral pressure (xy plane) and the system is evacuated in the z direction, the only thing that you might expect is that your system would be squeezed in that direction, and then the lateral pressure would relax. If you're thinking about a gas exerting a pressure on the interface, your system would need to have a pressurized gas instead of vacuum in the z direction. Please note that pressure would arise from the gas you have put into the box and not from an external pressure bath (I would stick to NVT for such a model). Among the possible issue you might face, the gas would probably be partially miscible in the liquid phase for high pressures. I hope it helps. cheers Andre On Wed, Apr 3, 2013 at 11:07 PM, Elisabeth katesed...@gmail.com wrote: Hello all, Does anyone know how one can study the effect of pressure on surface tension of pure liquids? Thanks, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- _ Prof. Dr. André Farias de Moura Department of Chemistry Federal University of São Carlos São Carlos - Brazil phone: +55-16-3351-8090 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- _ Prof. Dr. André Farias de Moura Department of Chemistry Federal University of São Carlos São Carlos - Brazil phone: +55-16-3351-8090 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: surface tension-crash at step 0
Hi Vitaly, To get surface tension one needs to extend Z to get a surface. From bulk it is not possible to create a surface. My problem is NVT crashes at step 0 and no energies or any other output file is written... On 5 April 2013 10:25, Dr. Vitaly Chaban vvcha...@gmail.com wrote: Hi all, I have problem ruunning NVT with Z direction extended to get surface tension of a polymer packed in a cell. I tried the same procedure for an alkane and NVT works well however for the polymer the run crashes at the very first step.. Doe anyone have clue what wrong could be? Many thanks. starting mdrun 'Polymer' 1000 steps, 1.0 ps. step 0 [node6:14338] *** Process received signal *** Hmmm. Is it not easier to get surface tension from bulk simulation?? On your error, I believe you broke some bond when extending your box in Z direction. Check energies at step 0 to be completely sure. Dr. Vitaly Chaban -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: surface tension-crash at step 0
The structure is taken from the end of previous simulation which terminates successfully. The initial structure should not be the problem. What I do only edit the last line (Z direction) of this structure to create vacuum. This the output from log file. Charge group distribution at step 0: 81 188 161 159 158 57 72 84 Grid: 3 x 3 x 6 cells Initial temperature: 425.057 K Started mdrun on node 0 Fri Apr 5 17:16:37 2013 Step Time Lambda 00.00.0 Energies (kJ/mol) Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14 5.40484e+038.26652e+032.54277e+031.31894e+03 -6.08804e+02 LJ (SR) Coulomb (SR) PotentialKinetic En. Total Energy 8.31491e+051.40229e+038.49817e+054.04825e+094.04910e+09 Conserved En.Temperature Pressure (bar) 4.04910e+091.12122e+087.78980e+08 and here is the output.log from the previous run (bulk simulation) DD step 19998999 vol min/aver 0.918 load imb.: force 2.0% Step Time Lambda 10001.00.0 Energies (kJ/mol) Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14 5.39477e+038.19320e+032.63606e+031.32991e+03 -5.94134e+02 LJ (SR) Coulomb (SR) PotentialKinetic En. Total Energy -4.34831e+031.34243e+031.39539e+041.52157e+042.91697e+04 Temperature Pressure (bar) 4.21421e+025.21511e+02 DD step 1999 vol min/aver 0.866 load imb.: force 1.0% Step Time Lambda 20002.00.0 Writing checkpoint, step 2000 at Sat Mar 9 22:16:02 2013 Energies (kJ/mol) Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14 5.29484e+038.23267e+032.54182e+031.31784e+03 -6.08717e+02 LJ (SR) Coulomb (SR) PotentialKinetic En. Total Energy -4.35248e+031.37024e+031.37962e+041.55820e+042.93782e+04 Temperature Pressure (bar) 4.31565e+02 -9.33582e+02 == ### == A V E R A G E S == ### == Statistics over 2001 steps using 201 frames Energies (kJ/mol) Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14 5.25008e+038.19447e+032.71542e+031.36750e+03 -5.69335e+02 LJ (SR) Coulomb (SR) PotentialKinetic En. Total Energy -4.29465e+031.31685e+031.39803e+041.53404e+042.93207e+04 Temperature Pressure (bar) 4.24873e+021.57552e-01 On 5 April 2013 16:34, Dr. Vitaly Chaban vvcha...@gmail.com wrote: But you have conf.gro On Fri, Apr 5, 2013 at 10:29 PM, Elisabeth katesed...@gmail.com wrote: Well I am not getting the trajectory. Simulation crashes before any step is calculated :( On 5 April 2013 11:20, Dr. Vitaly Chaban vvcha...@gmail.com wrote: My problem is NVT crashes at step 0 and no energies or any other output file is written... What about visualizing your interface? On 5 April 2013 10:25, Dr. Vitaly Chaban vvcha...@gmail.com wrote: Hi all, I have problem ruunning NVT with Z direction extended to get surface tension of a polymer packed in a cell. I tried the same procedure for an alkane and NVT works well however for the polymer the run crashes at the very first step.. Doe anyone have clue what wrong could be? Many thanks. starting mdrun 'Polymer' 1000 steps, 1.0 ps. step 0 [node6:14338] *** Process received signal *** Hmmm. Is it not easier to get surface tension from bulk simulation?? On your error, I believe you broke some bond when extending your box in Z direction. Check energies at step 0 to be completely sure. Dr. Vitaly Chaban -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] surface tension
Dear Dr. Moura, Thank you for your answer. I equilibrated the cell under NPT and then extended the Z direction to get a surface. Since I run NVT and Z is extended, I see no pressure dependence and surface tension values are similar for all P from 50 to 1000 bar. which is attributed to the fact that Z is extended and system tends to expand regardless of the pressure that was imposed in the initial NPT runs. I have no gas, it is liqiuid-vaccum and wanted to know if there is anyway to capture effect of pressure. From your answer I think there is no such possibility... Please comment if the answer is positive, Thank you On 4 April 2013 10:32, André Farias de Moura mo...@ufscar.br wrote: There's no simple answer for that. If you apply a lateral pressure (xy plane) and the system is evacuated in the z direction, the only thing that you might expect is that your system would be squeezed in that direction, and then the lateral pressure would relax. If you're thinking about a gas exerting a pressure on the interface, your system would need to have a pressurized gas instead of vacuum in the z direction. Please note that pressure would arise from the gas you have put into the box and not from an external pressure bath (I would stick to NVT for such a model). Among the possible issue you might face, the gas would probably be partially miscible in the liquid phase for high pressures. I hope it helps. cheers Andre On Wed, Apr 3, 2013 at 11:07 PM, Elisabeth katesed...@gmail.com wrote: Hello all, Does anyone know how one can study the effect of pressure on surface tension of pure liquids? Thanks, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- _ Prof. Dr. André Farias de Moura Department of Chemistry Federal University of São Carlos São Carlos - Brazil phone: +55-16-3351-8090 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] surface tension
Hello all, Does anyone know how one can study the effect of pressure on surface tension of pure liquids? Thanks, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] g_density: Segmentation fault
Dear all, I am trying to get the density profile for my liquid-vacuum interface using g_density -f trr -s tpr however g_density gives Segmentation fault. Does anyone had clue what could be wrong? Please comment, Thanks. Group 0 ( System) has XXX elements Group 1 ( Other) has XXX elements Select a group: 0 Selected 0: 'System' trn version: GMX_trn_file (single precision) Last frame 10 time 1.000 Segmentation fault -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] g_density: Segmentation fault
Hi Justin, Do I have to read less frames to circumvent the problem? I know g_density has been used for this purpose so there should be a way to resolve this. I am reading 1000 frames... Thanks On 2 April 2013 12:40, Justin Lemkul jalem...@vt.edu wrote: On 4/2/13 11:57 AM, Elisabeth wrote: Dear all, I am trying to get the density profile for my liquid-vacuum interface using g_density -f trr -s tpr however g_density gives Segmentation fault. Does anyone had clue what could be wrong? Please comment, Thanks. Group 0 ( System) has XXX elements Group 1 ( Other) has XXX elements Select a group: 0 Selected 0: 'System' trn version: GMX_trn_file (single precision) Last frame 10 time 1.000 Segmentation fault Segmentation faults are generic memory errors. Without a debugging backtrace, there's very little to suggest. -Justin -- ==**== Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justinhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ==**== -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/** Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/**Support/Mailing_Listshttp://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: density profile
Hi vitaly, The initial structure is indeed extended but the final output.gro is not. I think its because I am using the cpt file from the previous NPT runs as input for the new runs? Do I have to remove the -t flag? On 1 April 2013 12:47, Dr. Vitaly Chaban vvcha...@gmail.com wrote: Hi Elisabeth - The only explanation is that you actually DID NOT extend the box in Z direction. Look at the last line of confout.gro. g_density -d Z gives you a [local] density versus Z coordinate. Dr. Vitaly Chaban On Mon, Apr 1, 2013 at 5:33 PM, Elisabeth katesed...@gmail.com wrote: Hi Vitaly, I did NVT simulations and tried to obtain density profile at interface along Z using g_density -f .trr -s .tpr -d Z but I what I see is the density profile in the box not the interface. Box size is 3 nm and Before NVT runsI extended Z to 6 nm. Please see the attached profile. Thanks! -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: density profile
Hi Vitaly, The problem was with cpt file since it re sets the last line of gro. I removed the -f flag and now the Z direction is extended. However, I see that molecules tend to fill up the upper zone (free space) rapidly. I am wondering how I can obtain the density profile if I am going to get another uniformly distributed box after this NVT run? I am expecting to see how density changes with Z at the solvent -vacuum interface Please advise me on this,, Thanks! On 1 April 2013 13:14, Dr. Vitaly Chaban vvcha...@gmail.com wrote: I think if you use checkpoint files, the program does not read either MDP, or GRO, or TOP, or anything except CPT. Dr. Vitaly Chaban On Mon, Apr 1, 2013 at 7:10 PM, Elisabeth katesed...@gmail.com wrote: Hi vitaly, The initial structure is indeed extended but the final output.gro is not. I think its because I am using the cpt file from the previous NPT runs as input for the new runs? Do I have to remove the -t flag? On 1 April 2013 12:47, Dr. Vitaly Chaban vvcha...@gmail.com wrote: Hi Elisabeth - The only explanation is that you actually DID NOT extend the box in Z direction. Look at the last line of confout.gro. g_density -d Z gives you a [local] density versus Z coordinate. Dr. Vitaly Chaban On Mon, Apr 1, 2013 at 5:33 PM, Elisabeth katesed...@gmail.com wrote: Hi Vitaly, I did NVT simulations and tried to obtain density profile at interface along Z using g_density -f .trr -s .tpr -d Z but I what I see is the density profile in the box not the interface. Box size is 3 nm and Before NVT runsI extended Z to 6 nm. Please see the attached profile. Thanks! -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] density profile
Dear all, In order to get density profile of a pure system the box has to get extended in one drection (e.g Z) as we do for surface tension calculations or density profile can be also obtained from the usual box filled up with the molecules without need for empty space in Z Thanks for any comments and advise in advance :) Best E. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] density profile
Thank Justin. I am interested in the density profile at the solvent-air interface. I wanted to see how density changes with position at different pressures...I have the equilibrated boxes at several pressures obtained from NPT but I am not sure whether running g_density on the current simulations cells provides the profile I am after. Do you have any clue whether or not the cell has to be extended in one direction to obtain the density profile? Thanks! On 31 March 2013 12:21, Justin Lemkul jalem...@vt.edu wrote: On 3/31/13 11:41 AM, Elisabeth wrote: Dear all, In order to get density profile of a pure system the box has to get extended in one drection (e.g Z) as we do for surface tension calculations or density profile can be also obtained from the usual box filled up with the molecules without need for empty space in Z Thanks for any comments and advise in advance :) The density of the system is written to the .edr file. There is no need (or use) for a density profile if the system is homogeneous. -Justin -- ==**== Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justinhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ==**== -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/** Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/**Support/Mailing_Listshttp://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] density profile
Thanks Justin for your reply. To avoid compressing down the cell I thought using semiisotropic option with 0 compressibility in Z would be appropriate. Here are the steps: 1- First I did a 10 ns NPT to equilibrate the box and used the last frame gro and cpt file as input for the next step 2- I extended the box in Z to more than twice the initial box size and issued the following: grompp -f md.mdp -c extendedZ.gro -t .cpt -p .top -o extendedZ.tpr mpirun -np 4 mdrun_mpi -deffnm .extendedZ -s -o -c -g -e -x -v My questions are: 1) Can I get the proper solvent-air interface to obtain density profile by extending the Z direction (last line of the gro file obtained from NPT at each pressure)? Is that what yo mean by then build a new unit cell and run under NVT ? 2) Is the semi isotropic option with 0 compressibility in Z is appropriate to keep the pressure fixed? (as shown below in mdp settings)? 3) If the answer to q 2 is negative, whats is appropriate way to shift from NPt to NVT run. Do I have to use the equilibrated density from NPT runs and edit the last line of gro from -c flag of NPT.gro to obtain the correct density? The reason I am asking is that since P fluctuations are huge I am not sure if the the density of box from last frame of -c NPT.gro output is the correct density... Please advise on the above queries..Thanks a lot Below is the md.mdp contents: pbc = xyz integrator = md dt = 0.001 nsteps = 2000 nstcomm = 100 ; Output control nstenergy = 100 nstxout = 100;1 nstvout = 0 nstfout = 0 nstlog = 1000 nstxtcout = 1000 ; Neighbor searching nstlist = 10 ns_type = grid ; Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 1 rvdw-switch = 1 ;0.9 ; Cut-offs rlist = 1.35 rcoulomb= 1.1 ;1.1 rvdw= 1.1 Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 425 * * *Pcoupl = berendsen* *Pcoupltype = semiisotropic * *tau_p = 1 0.5* *compressibility = 3.5e-5 0* *ref_p = 50 ** * gen_vel = no gen_temp= 500.0 gen_seed= 173529 constraints = none constraint-algorithm = lincs On 31 March 2013 12:47, Justin Lemkul jalem...@vt.edu wrote: On 3/31/13 12:27 PM, Elisabeth wrote: Thank Justin. I am interested in the density profile at the solvent-air interface. I wanted to see how density changes with position at different pressures...I have the equilibrated boxes at several pressures obtained from NPT but I am not sure whether running g_density on the current simulations cells provides the profile I am after. Do you have any clue whether or not the cell has to be extended in one direction to obtain the density profile? If you're trying to produce an air-water interface, then indeed you will need some model for air within the unit cell, but there are several practical problems with that. The first is that g_density does not deal well with changing box vectors, and the density profiles it produces under an NPT ensemble are incorrect (outstanding bug that needs to be fixed, IIRC). The second is that if you run a simulation under NPT with void space, the unit cell will just compress down until there is no empty space. Probably the best solution would be to equilibrate under NPT at the desired conditions, then build a new unit cell and run under NVT. Using NVT will circumvent both problems, and if the proper pressure, density, etc have already been achieved, then you've got what you're after. -Justin -- ==**== Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justinhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ==**== -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/** Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/**Support/Mailing_Listshttp://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive
Re: [gmx-users] Re: density profile
Thanks Vitaly. I am wondering what is the use of semiisotropic with 0 compressibility in Z then? I was hoping to run NPT to secure a fixed pressure. I also wanted to know if surface tension can be also calculated under NVT (if NPT fails for this puporse) Thanks! On 31 March 2013 14:01, Dr. Vitaly Chaban vvcha...@gmail.com wrote: Thanks Justin for your reply. To avoid compressing down the cell I thought using semiisotropic option with 0 compressibility in Z would be appropriate. You must use NVT only. Otherwise, the cell will compress in XY directions to compensate its inability to compress in Z direction. Here are the steps: 1- First I did a 10 ns NPT to equilibrate the box and used the last frame gro and cpt file as input for the next step 2- I extended the box in Z to more than twice the initial box size and issued the following: grompp -f md.mdp -c extendedZ.gro -t .cpt -p .top -o extendedZ.tpr mpirun -np 4 mdrun_mpi -deffnm .extendedZ -s -o -c -g -e -x -v My questions are: 1) Can I get the proper solvent-air interface to obtain density profile by extending the Z direction (last line of the gro file obtained from NPT at each pressure)? Is that what yo mean by then build a new unit cell and run under NVT ? If you want solvent/air interface, you should add $air. If you want solvent/vacuum interface, it is enough to extend the box. If you want liquid/vapor interface, you need to re-equilibrate the system in NVT with a space for vapor available in the box. Dr. Vitaly Chaban -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] semi isotropic settings
Hi all, I wanted to use Pcoupltype = semiisotropicbut I am not clear about the compressibility settings if I want to keep Z directions fixed in my system and letting X and Y change. Which case 1 or 2 is correct for semi isotropic option? 1) compressibility = 3.5e-5 3.5e-5 0 2 ) or compressibility = 1 0? Manual says 2 values are needed for x/y and z directions respectively. Thanks -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] semi isotropic settings
Hi Justin, The corect value is 3.5e-5 but I thoiugh x/y means the ratio od these two. If thats nto the case I guess compressibility = 3.5e-5 0 is the correct value. What do you think? On 30 March 2013 11:41, Justin Lemkul jalem...@vt.edu wrote: On 3/30/13 11:39 AM, Elisabeth wrote: Hi all, I wanted to use Pcoupltype = semiisotropicbut I am not clear about the compressibility settings if I want to keep Z directions fixed in my system and letting X and Y change. Which case 1 or 2 is correct for semi isotropic option? 1) compressibility = 3.5e-5 3.5e-5 0 2 ) or compressibility = 1 0? Manual says 2 values are needed for x/y and z directions respectively. Case #2, though I doubt a value of 1 is a correct compressibility. But as the manual says, you need two values - one for x/y (i.e. those two dimensions are equally compressible) and one for z. -Justin -- ==**== Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justinhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ==**== -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/** Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/**Support/Mailing_Listshttp://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] surface tension-density profile
Deal all, I am trying to build up my alkane system to calculate the surface tension a fixed pressure. Here are the steps: 1- First I did a 10 ns NPT to equilibrate the box and used the last frame gro and cpt file as input for the next step 2- I extended the box in Z to more than twice the initial box size and issued the following: grompp -f md.mdp -c extendedZ.gro -t .cpt -p .top -o extendedZ.tpr mpirun -np 4 mdrun_mpi -deffnm .extendedZ -s -o -c -g -e -x -v However, I am not clear about the setting of pressure coupling and compressibility. I guess I need to use 0 compressibility in Z with isotropic option? Please advise me on the settings and details of these calculations. I am going to compare the surface tension and density profiles at several different pressures. Thank you in advance :) Below is the md.mdp contents: pbc = xyz integrator = md dt = 0.001 nsteps = 2000 nstcomm = 100 ; Output control nstenergy = 100 nstxout = 100;1 nstvout = 0 nstfout = 0 nstlog = 1000 nstxtcout = 1000 ; Neighbor searching nstlist = 10 ns_type = grid ; Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 1 rvdw-switch = 1 ;0.9 ; Cut-offs rlist = 1.35 rcoulomb= 1.1 ;1.1 rvdw= 1.1 Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 425 * * *Pcoupl = berendsen* *Pcoupltype = isotropic * *tau_p = 1 0.5* *compressibility = 3.5e-5 3.5e-5 0* *ref_p = 50 10** * gen_vel = no gen_temp= 500.0 gen_seed= 173529 constraints = none constraint-algorithm = lincs -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Shift functions
On 10 February 2012 09:41, Mark Abraham mark.abra...@anu.edu.au wrote: On 11/02/2012 1:19 AM, Elisabeth wrote: Hello all, Does the shift function use group based truncation? See the discussion of charge groups in manual section 3.4.2. Thanks Mark. -1- First of all if I am right charge groups in gromacs language in identical to group based truncations? Using charge groups as the indivisible entity upon which neighbour list construction is based is using group based truncations. The usual alternative is using atoms as the, well, *atomic* unit. The former can be equivalent to the latter if there's one atom per charge group. Manual 342: This reduces the cut-off effects from the charge-charge level to the dipole-dipole level, which decay much faster 2- I am not able to realize why we go from charge-charge the dipole-dipole changes? Charge groups are constructed to have neutral charge (or integer charge where necessary). To first order, the difference between any such group being in the neighbour list of another such group or not (according to the cut-off radius) is equivalent to a point dipole being in the cut-off sphere or not. Charge groups with arbitrary charge (or partially-charged atoms, when using atom-based truncation) do not have this quality, and the distance at which a charge-charge interaction is significant is much larger than that of a dipole-dipole. In the manual I see: by using shifted forces there is no need for charge groups (=group based?!) in the neighbor list? Can anyone shed some light on calculation of shifted forces? What's not clear from the above and 4.1.5? 3- I understand that use of shift makes the potentials have continuous derivatives at cutoffs but that how this makes use of charge groups unnecessary, I dont see! Remember that the neighbour list is constructed to permit the computation of a finite number of interactions with the central atom/group. You want to stop computing them when they're close enough to zero that you don't care. If they actually go to zero, then you don't care. If they decay as 1/r (charge-charge) then at typical r_c values you should care. If they decay as 1/r^2 (dipole-dipole, IIRC) then at typical r_c values things are OK. If the value of the force is non-zero at the cut-off, then there is an interaction at that distance and not one just past that distance. This generates artefacts that are serious for non-zero charges at the kinds of cut-off distances for which force fields are parametrized, but much less serious if computed over neutral charge groups. If the value of the force is zero at the cut-off (i.e. shift potential), then no atom or charge group has any interaction with the central group/atom at that distance, so you don't need to care about whether the truncation is based on atoms or groups. You do have to care about the effect of the modified Coulomb potential, however. Hello Mark, I read over your answer several times. I am still unclear about modified Coulomb potential. In the manual modified Coulomb potential refers to shift/switch functions, as I realized. so when the force is zero at r_c and it doesnt matter whether truncation is based on atoms or groups, why effect of modified coulomb potential is important? -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Shift functions
Hello all, Does the shift function use group based truncation? In the manual I see: by using shifted forces there is no need for charge groups (=group based?!) in the neighbor list? Can anyone shed some light on calculation of shifted forces? Thanks, Eli -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Shift functions
Hello all, Does the shift function use group based truncation? See the discussion of charge groups in manual section 3.4.2. Thanks Mark. -1- First of all if I am right charge groups in gromacs language in identical to group based truncations? Manual 342: This reduces the cut-off effects from the charge-charge level to the dipole-dipole level, which decay much faster 2- I am not able to realize why we go from charge-charge the dipole-dipole changes? In the manual I see: by using shifted forces there is no need for charge groups (=group based?!) in the neighbor list? Can anyone shed some light on calculation of shifted forces? What's not clear from the above and 4.1.5? 3- I understand that use of shift makes the potentials have continuous derivatives at cutoffs but that how this makes use of charge groups unnecessary, I dont see! 4- and based on 3, shift forces dont neglect tail corrections for LJ as cutoffs do? Am I correct? Thank you -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Equilibration of polymer chains- gauche defects
Hi all, Can you please guide me how one can calculate fraction of gauche+ and gauche- defects in the center of a polymer chain? Appreciate you help. Juliete -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] phase transition in MD simulation
Dear experts, I am working on a pure polymer melt system under NPT and at a given T, pressure is increased up to around 1000 bar. However, the phase diagram of the polymer is indicating that at this T, as pressure is increased to above 500 bar, system falls below the melting point. I mean I am doing NPT on a liquid polymer by packing the chains in to the box but even applying a pressure of 1000 bar does not make the atoms or chains resemble a solid phase. The system would only become a more packed liquid phase. In other words, MD is representing an extrapolated liquid phase in the subcooled region. My question is whether there is any way to identify phase transition from MD simulation. Please comment and let me know what you think. Thanking you, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] index file subgroups
Hello, I am trying to figure out how to plot or obtain the average values of Gyration radius from different subgroups of index file. index file prompts for a single group and I need to calculate the averages for all subgroup by hand. Is there any way to compute averages or plot all subgroups directly? g_polystat -f *.trr -s *.tpr -o *.xvg -n *.ndx -w Thanks, Regards, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] NPT - density off
On 20 August 2011 21:20, Mark Abraham mark.abra...@anu.edu.au wrote: On 20/08/2011 8:02 AM, Elisabeth wrote: Dear experts, I am intending to calculate the equilibrium density of a pure hydrocarbon at different pressures , at 425 K. The normal boiling point of the liquid is around 350 K. For pressures below 100 bar densities I am getting from NPT is in accurate. I start form a structure which is compressed to above 0.6 g/cm3 density but since temperature is high density goes than to the values below. P = 50 bar NPT rho= 0.344 experimental density ~ 0.54 P100 NPT rho= 0.43 experimental density ~ 0.55 densities become more accurate for P 100 bar. P 500 rho= 0.56experimental density ~ 0.61 I thought maybe you have some idea on how this inaccuracy can be improved. These could reflect limitations in the model you are using (which almost certainly wasn't parametrized upon data like this). Your integration time step is twice as large as is commonly used in the absence of constraints. Also, be sure you are measuring your density only after equilibration, not an average that includes the equilibration period - and describe that method so people here know you're making such sensible measurements. Hello Mark, Thank you for your comment. Densities are equilibrated. I am not sure what you mean by limitations in the model. I actually fixed the density and tried NVT using the experimental density at 298 and found a full agreement between vaporization heat of the liquid and that reported in the original OPLS paper. Once I try to obtain the density using NPT this inaccuracy in densities appear. I have been using 1 fs dt with berendsen barostat for 5 ns and collected data over the last 1 ns ( there are about 2500 atoms in the system). So now I am going to follow your instruction and try dt of 0.5 fs for another 5ns along with cpt files from previous barostat trial (berendsen). I will get back to you to report the results soon. BTW: I am using version 4.5.4 Thank you, Mark Thank you in advance for your invaluable help. Best, ;Bonds constraints = none constraint-algorithm = lincs ;Run control integrator = md dt = 0.001 nsteps = 500 nstcomm = 100 ;Output control nstenergy = 100 nstxout = 100 nstvout = 0 nstfout = 0 nstlog = 1000 nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 425 ;Pressure coupling Pcoupl = berendsen Pcoupltype = isotropic tau_p = 1 compressibility = 3.5e-5 ref_p = 100 ;Velocity generation gen_vel = no;yes gen_temp= 425 gen_seed= 173529 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] NPT - density off
Dear experts, I am intending to calculate the equilibrium density of a pure hydrocarbon at different pressures , at 425 K. The normal boiling point of the liquid is around 350 K. For pressures below 100 bar densities I am getting from NPT is in accurate. I start form a structure which is compressed to above 0.6 g/cm3 density but since temperature is high density goes than to the values below. P = 50 bar NPT rho= 0.344 experimental density ~ 0.54 P100 NPT rho= 0.43 experimental density ~ 0.55 densities become more accurate for P 100 bar. P 500 rho= 0.56experimental density ~ 0.61 I thought maybe you have some idea on how this inaccuracy can be improved. Thank you in advance for your invaluable help. Best, ;Bonds constraints = none constraint-algorithm = lincs ;Run control integrator = md dt = 0.001 nsteps = 500 nstcomm = 100 ;Output control nstenergy = 100 nstxout = 100 nstvout = 0 nstfout = 0 nstlog = 1000 nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 425 ;Pressure coupling Pcoupl = berendsen Pcoupltype = isotropic tau_p = 1 compressibility = 3.5e-5 ref_p = 100 ;Velocity generation gen_vel = no;yes gen_temp= 425 gen_seed= 173529 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] berendsen P coupling and fluctuation properties
Hello, Thank you. I am looking at potential energies to calculate vaporization heat. I wanted to know how the fact that berendsen does not lead to correct ensemble is affecting total potential energy of the system. I have an unclear image of correct ensemble. Does this mean whatever output I am getting from my runs are unreliable? Appreciate any clarification. Thank you, Best, On 15 August 2011 19:34, Mark Abraham mark.abra...@anu.edu.au wrote: On 16/08/2011 7:05 AM, Elisabeth wrote: Dear all, I noticed that applying Parrinello-Rahman (PR) pressure coupling even after equilibration with berendsen does not lead to target value for pressure when ;Bonds constraints = none is used. The use of constraints and the integration step size is linked. Roughly speaking, no constraints should accompany a 0.5 fs time step, H-bond constraints with 1fs and all constraints with 2fs. Haphazard changes to .mdp files have all kinds of these gotchas. I tried Berendsen for to get fixed pressure ( 50 bar) but in the next run PR even for long time is giving 52 bar. This is the case for other target pressures too. You need to be sure to collect statistics only after equilibration, and consider whether the observed variation is consistent with convergence to a given value. So this made me select Berenden which is giving target pressure values but my concern is whether my results are reliable because BR does not give the exact ensemble as PR. I read somewhere on the list that fluctuation properties can not be calculated when BR is used. What does fluctuation property mean? BR does not produce the correct ensemble. I forget the details about why, but there are references in the T-coupling section of the manual you should consider. Does this mean that any property calculated form fluctuations of some other quantity can not be obtained? like heat capacity which is defined based on enthalpy fluctuations? IIRC, yes. I am interested in potential energy terms (g_energy bonded/non boned terms) I routinely struggle to see why people think they can learn anything from these. and structural properties like rdf for a number of polymer molecules, system size around 3000 atoms. Thank you for your comments. Best, constraints = none ;Run control integrator = md dt = 0.001 nsteps = 500 nstcomm = 100 nstenergy = 100 nstxout = 100 nstlist = 10 ns_type = grid coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 Pcoupl = berendsen Pcoupltype = isotropic tau_p = 1 compressibility = 3.5e-5 ref_p = 10 gen_vel = no gen_temp= 300.0 gen_seed= 173529 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] berendsen P coupling and fluctuation properties
Thanks Justin for your thorough description. Is there any specific tool in gromacs to plot histograms or you mean I follow the customary plotting procedure? I have one more question. I am now trying using PR pressure scheme with v-rescale thermostat for NPT dynamics. Before dynamics run, can I apply berendsen barostat combined with v-rescale thermostat to equilibrate or since I am using berendsen barostat, T coupling scheme has to be also berendsen? Thanks, On 16 August 2011 13:17, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote:f Hello, Thank you. I am looking at potential energies to calculate vaporization heat. I wanted to know how the fact that berendsen does not lead to correct ensemble is affecting total potential energy of the system. I have an unclear image of correct ensemble. Does this mean whatever output I am getting from my runs are unreliable? The Berendsen coupling algorithms produce very narrow distributions of temperature and pressure. These do not correspond to the correct distribution of a true statistical mechanical ensemble; plot a histogram and you'll find that the results are shockingly different between Berendsen and, say, Nose-Hoover for T. Therefore, the ensemble you're using to measure your properties is not NPT (or NVT, in the case of thermostat only), it's something undefined and likely not real. Applying better thermostats and barostats that have been shown to produce the correct distributions is more rigorously correct. The effects are usually not terribly noticeable unless you go looking for them. You'll find that the Berendsen algorithms quite faithfully give you the target T and P that you desire (which does make them very good for initial equilibration) and so you naturally assume that everything is fine. While that's all well and good in one sense, the fluctuations of T (and thus velocities, and thus kinetic energy) are wrong. If KE does not fluctuate properly, then neither does PE (since you're then affecting energy flow between PE and KE, while total energy stays fixed, in theory). The larger point is that if you claim to apply an NVT (or NPT) ensemble using Berendsen coupling, you're not correct, strictly speaking. -Justin Appreciate any clarification. Thank you, Best, On 15 August 2011 19:34, Mark Abraham mark.abra...@anu.edu.au mailto: mark.abra...@anu.edu.**au mark.abra...@anu.edu.au wrote: On 16/08/2011 7:05 AM, Elisabeth wrote: Dear all, I noticed that applying Parrinello-Rahman (PR) pressure coupling even after equilibration with berendsen does not lead to target value for pressure when ;Bonds constraints = none is used. The use of constraints and the integration step size is linked. Roughly speaking, no constraints should accompany a 0.5 fs time step, H-bond constraints with 1fs and all constraints with 2fs. Haphazard changes to .mdp files have all kinds of these gotchas. I tried Berendsen for to get fixed pressure ( 50 bar) but in the next run PR even for long time is giving 52 bar. This is the case for other target pressures too. You need to be sure to collect statistics only after equilibration, and consider whether the observed variation is consistent with convergence to a given value. So this made me select Berenden which is giving target pressure values but my concern is whether my results are reliable because BR does not give the exact ensemble as PR. I read somewhere on the list that fluctuation properties can not be calculated when BR is used. What does fluctuation property mean? BR does not produce the correct ensemble. I forget the details about why, but there are references in the T-coupling section of the manual you should consider. Does this mean that any property calculated form fluctuations of some other quantity can not be obtained? like heat capacity which is defined based on enthalpy fluctuations? IIRC, yes. I am interested in potential energy terms (g_energy bonded/non boned terms) I routinely struggle to see why people think they can learn anything from these. and structural properties like rdf for a number of polymer molecules, system size around 3000 atoms. Thank you for your comments. Best, constraints = none ;Run control integrator = md dt = 0.001 nsteps = 500 nstcomm = 100 nstenergy = 100nstxout = 100 nstlist = 10 ns_type = grid coulombtype = Shiftvdw-type = Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist
[gmx-users] berendsen P coupling and fluctuation properties
Dear all, I noticed that applying Parrinello-Rahman (PR) pressure coupling even after equilibration with berendsen does not lead to target value for pressure when ;Bonds constraints = none is used. I tried Berendsen for to get fixed pressure ( 50 bar) but in the next run PR even for long time is giving 52 bar. This is the case for other target pressures too. So this made me select Berenden which is giving target pressure values but my concern is whether my results are reliable because BR does not give the exact ensemble as PR. I read somewhere on the list that fluctuation properties can not be calculated when BR is used. What does fluctuation property mean? Does this mean that any property calculated form fluctuations of some other quantity can not be obtained? like heat capacity which is defined based on enthalpy fluctuations? I am interested in potential energy terms (g_energy bonded/non boned terms) and structural properties like rdf for a number of polymer molecules, system size around 3000 atoms. Thank you for your comments. Best, constraints = none ;Run control integrator = md dt = 0.002 nsteps = 500 nstcomm = 100 nstenergy = 100 nstxout = 100 nstlist = 10 ns_type = grid coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 Pcoupl = berendsen Pcoupltype = isotropic tau_p = 1 compressibility = 3.5e-5 ref_p = 10 gen_vel = no gen_temp= 300.0 gen_seed= 173529 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: Heat of Vaporization
On 3 August 2011 15:31, Dr. Vitaly V. Chaban vvcha...@gmail.com wrote: On Wed, Aug 3, 2011 at 3:17 PM, Elisabeth katesed...@gmail.com wrote: On 2 August 2011 15:29, Dr. Vitaly V. Chaban vvcha...@gmail.com wrote: Hello, I wanted to know your ideas on calculation of heat of vaporization using a single phase run rather than running two separate simulations for liquid and gas! 1- Two separate simulations for liquid and gas DHvap = Ugas - Uliq + RT 1a: total liquid potential - *total* potential of a single chain in vacu ( bond+angle+torsion + nonbonded interaction of chain with itself) or 1b: *total liquid potential - intra potential of a single chain in vacu ( bond+angle+torsion) * 2- Single liquid phase run: (non need to run in vacu) 2a : DHvap = total liquid potential - intra molecular potential terms in liquid phase (same liquid phase simulation by adding up bond+angle+torsion terms) 2a: In other worlds *DHvap= Uliq-nonbonded (vdw+electrostatics) * In my case the latter definition is giving much more accurate results than 1a. I would like to know your idea and comments on methods 1b and 2a. Appreciate your comments. What is your system? system is liquid hydrocarbon polymer and 1a is giving inaccurate values. 2a works far better but it seems not to be a common method. First of all, you should understand which particles exist in the vapor phase of your polymer. Notwithstanding the atomistic simulation. If this question is answered correctly, any method will provide you a decent result. I don't think that the whole chains of this polymer are flying in the vapor phase. Hello, Thanks. I have only one single chain in vauo. pbc = no and cutoff is set to zero for method 1a. Simulation runs for 20 ns and is equilibrated. Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Heat of Vaporization
Hello, I wanted to know your ideas on calculation of heat of vaporization using a single phase run rather than running two separate simulations for liquid and gas! 1- Two separate simulations for liquid and gas DHvap = Ugas - Uliq + RT 1a: total liquid potential - *total* potential of a single chain in vacu ( bond+angle+torsion + nonbonded interaction of chain with itself) or 1b: *total liquid potential - intra potential of a single chain in vacu ( bond+angle+torsion) * 2- Single liquid phase run: (non need to run in vacu) 2a : DHvap = total liquid potential - intra molecular potential terms in liquid phase (same liquid phase simulation by adding up bond+angle+torsion terms) 2a: In other worlds *DHvap= Uliq-nonbonded (vdw+electrostatics) * In my case the latter definition is giving much more accurate results than 1a. I would like to know your idea and comments on methods 1b and 2a. Appreciate your comments. Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] OPLSAA parameters
Hi Justin, So far I have been #including ffoplsaa.itp and running simulations. The reason I was after parameters is just to have a better indea on the FF parameters rather than just using them blindly, that's it :) Thank you, On 31 July 2011 20:03, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: Thanks Justin for your guidance. By following your tips, I have now all the parameters I need... Please let me know if I am missing any [ angletypes ] or [ dihedraltypes ] for a hydrocarbon molecule. Seems fine to me; grompp will complain if something is missing. Is there a reason that you're trying to parse out certain parameters? It'd be far easier to just #include oplsaa.ff/forcefield.itp in your topology and have the entire OPLS-AA force field available to you, rather than reinventing the wheel. -Justin Best, [ bondtypes ] ; i j func b0 kb CTCT 10.15290 224262.4 ; CHARMM 22 parameter file CTHC 10.10900 284512.0 ; CHARMM 22 parameter file [ angletypes ] ;i j k func th0 cth CT CT CT 1 112.700488.273 ; CHARMM 22 parameter file CT CT HC 1 110.700313.800 ; CHARMM 22 parameter file HC CT HC 1 107.800276.144 ; CHARMM 22 parameter file [ dihedraltypes ] CT CT CT CT 3 2.92880 -1.46440 0.20920 -1.67360 0.0 0.0 ; hydrocarbon all-atom CT CT CT HC 3 0.62760 1.88280 0.0 -2.51040 0.0 0.0 ; hydrocarbon all-atom HC CT CT HC 3 0.62760 1.88280 0.0 -2.51040 0.0 0.0 ; hydrocarbon *new* 11/99 On 30 July 2011 07:40, Justin A. Lemkul jalem...@vt.edu mailto: jalem...@vt.edu wrote: Elisabeth wrote: Hello to all, I am looking for bonded and nonbonded parameters for C H atoms in saturated hydrocarbons. atom types are: Carbon atoms opls_135 opls_136 H opls_140 opls_135 12.01100 ; alkane CH3 opls_136 12.01100 ; alkane CH2 opls_1401.00800 ; alkane H. Please help me out with the following: - I am wondering if I am selecting the correct values, for C-C and C-H bonds from ffoplsaabon.itp. The C-C bond is not correct for alkanes. It is for the C-C bond of oxalic acid, as the comment indicates. In ffnonbonded.itp, you can see how opls_* types are mapped to different bonded atom types (for instance, opls_135/136 are CT). Note that in the topology, you do not have to explicitly define values when using OPLS-AA; any bond between two defined bonded atom types will be automatically detected and used. C C 10.15100 292880.0 ; wlj oxalic acid, etc. HCC 10.10900 284512.0 ; wlj 7/96 - what does wlj mean? Probably someone's initials, which occur frequently to indicate either the source of the parameters or the individual that added them. In this case, my guess would be Jorgensen is the referenced individual here. - what are the units? This is in the manual. Angle CA CA CZ 1 120.000585.760 ; wlj CA CA CR 1 120.000527.184 ; CA CA CX 1 120.000711.280 ; C CT Cl 1 109.800577.392 ; wlj - I can not distinguish CA CZ What are they referring to? See the above tip about converting atom types in ffnonbonded.itp. [ dihedraltypes ] C C CT HC 3 0.17782 0.53346 0.0 -0.71128 0.0 0.0 ; dicarbonyls BMC 8,1881(2000) - Is this what I need for oplsaa 135 136 140? Probably not. The comment indicates the dihedral is for dicarbonyls. The appropriate dihedral is: CT CT CT HC 3 0.62760 1.88280 0.0 -2.51040 0.0 0.0 ; hydrocarbon all-atom Note how the comment indicates it is for use with hydrocarbons. Again, like bonds as described above, you do not need to explicitly list parameters in the topology. They will be inferred from the atom types. ffoplsaanb.itp. opls_135 CT6 12.01100-0.180 A 3.5e-01 2.76144e-01 opls_136 CT6 12.01100-0.120 A 3.5e-01 2.76144e-01 opls_140 HC1 1.00800 0.060 A 2.5e-01 1.25520e-01 - are these parameters for opls_135 12.01100 ; alkane CH3 opls_136 12.01100 ; alkane CH2 opls_1401.00800 ; alkane H. ? Yes. Also in usr/local/gromacs/share/__**gromacs/top/ directory I can not find ffoplsaanb,itp and ffoplsaabon.itp for 4.5.4 whereas in previous versions these files were accessible. I am
[gmx-users] citation
Hello all, Can anyone tell me how gromacs may be cited besides Principal Papers 1. Berendsen, et al. (1995) *Comp. Phys. Comm.* *91: *43-56. (DOIhttp://dx.doi.org/doi:10.1016/0010-4655%2895%2900042-E_, Citations of this paperhttp://scholar.google.com/scholar?cites=6703762643559347673hl=en ) 2. Lindahl, et al. (2001) *J. Mol. Model.* *7:* 306-317. (DOIhttp://dx.doi.org/10.1007/s008940100045, Citations of this paperhttp://scholar.google.com/scholar?cites=135270318718960hl=en ) 3. van der Spoel, et al. (2005) *J. Comput. Chem.* *26:* 1701-1718. (DOIhttp://dx.doi.org/10.1002/jcc.20291, Citations of this paperhttp://scholar.google.com/scholar?cites=8599454604037079587hl=en ) 4. Hess, et al. (2008) *J. Chem. Theory Comput.* *4:* 435-447. (DOIhttp://dx.doi.org/10.1021/ct700301q, Citations of this paperhttp://scholar.google.com/scholar?cites=15361472353040457616hl=en ) I am using oplsaa and found ; References for the OPLS-AA force field: ; ; W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, ; J. Am. Chem. Soc. 118, 11225-11236 (1996). ; W. L. Jorgensen and N. A. McDonald, Theochem 424, 145-155 (1998). ; W. L. Jorgensen and N. A. McDonald, J. Phys. Chem. B 102, 8049-8059 (1998). ; R. C. Rizzo and W. L. Jorgensen, J. Am. Chem. Soc. 121, 4827-4836 (1999). ; M. L. Price, D. Ostrovsky, and W. L. Jorgensen, J. Comp. Chem. (2001). ; E. K. Watkins and W. L. Jorgensen, J. Phys. Chem. A 105, 4118-4125 (2001). ; G. A. Kaminski, R.A. Friesner, J.Tirado-Rives and W.L. Jorgensen, J. Phys. Chem. B 105, 6474 (2001). ; Do I have to include all these or some of them can be skipped? Citations depend on specific features ( e.g free energy calculations, ...) that have been utilized or the above references are sufficient? Am I missing anything? Thank you, Best -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] OPLSAA parameters
Hello to all, I am looking for bonded and nonbonded parameters for C H atoms in saturated hydrocarbons. atom types are: Carbon atoms opls_135 opls_136 H opls_140 opls_135 12.01100 ; alkane CH3 opls_136 12.01100 ; alkane CH2 opls_1401.00800 ; alkane H. Please help me out with the following: - I am wondering if I am selecting the correct values, for C-C and C-H bonds from ffoplsaabon.itp. C C 10.15100 292880.0 ; wlj oxalic acid, etc. HCC 10.10900 284512.0 ; wlj 7/96 - what does wlj mean? - what are the units? Angle CA CA CZ 1 120.000585.760 ; wlj CA CA CR 1 120.000527.184 ; CA CA CX 1 120.000711.280 ; C CT Cl 1 109.800577.392 ; wlj - I can not distinguish CA CZ What are they referring to? [ dihedraltypes ] C C CT HC 3 0.17782 0.53346 0.0 -0.71128 0.0 0.0 ; dicarbonyls BMC 8,1881(2000) - Is this what I need for oplsaa 135 136 140? ffoplsaanb.itp. opls_135 CT6 12.01100-0.180 A3.5e-01 2.76144e-01 opls_136 CT6 12.01100-0.120 A3.5e-01 2.76144e-01 opls_140 HC1 1.00800 0.060 A2.5e-01 1.25520e-01 - are these parameters for opls_135 12.01100 ; alkane CH3 opls_136 12.01100 ; alkane CH2 opls_1401.00800 ; alkane H. ? Also in usr/local/gromacs/share/gromacs/top/ directory I can not find ffoplsaanb,itp and ffoplsaabon.itp for 4.5.4 whereas in previous versions these files were accessible. I am using itp files from version 4.0.7 and wanted to make sure the above parameters for oplsa_135 oplsa_136 and oplsa_140 have been remained unchanged in 4.5.4. Where can I find these files? Thanks, Regards, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] error: Only triclinic boxes...
On 3 June 2011 03:25, Mark Abraham mark.abra...@anu.edu.au wrote:
On 3/06/2011 3:06 AM, Elisabeth wrote:
Hello all,
I am getting the error below at the very beginning of the simulation (both
serial and parallel). I am sure I did not encounter this problem before with
the same input files. This has just happened now. I really have no clue why
this is happening. could you please help me? Thank you all in advance.
You broke the simulation. Do not use parinello-rahman for equilibration.
It's only good close to equilibrium, and can oscillate wildly under the
wrong conditions. Generate velocities on a structure whose density is
appropriate for the ensemble you're working with, equilibrate, then compress
in stages, very gently. I've lost track of the number of times I've
suggested people not just randomly apply some large pressure on some random
configuration and expect it to work.
Dear gmx experts,
The reason I was getting error below
Warning: Only triclinic boxes with the first vector parallel to the x-axis
and the second vector in the xy-plane are supported.
Box (3x3):
Box[0]={ nan, 0.0e+00, 0.0e+00}
Box[1]={ nan, nan, nan}
Box[2]={ nan, nan, nan}
Can not fix pbc.
was that I had mistakenly decreased box size to a size considerably smaller
than initial structure. One of my colleagues provided me with input files to
get an idea how these runs work and we were surprised how come the
simulation that would run before I receive it, is giving error message
above. I checked on gro file several times as you suggested but did not
notice the last line..!
Thanks Mark, Justin and Dr. Chaban for your comments.
Best,
Warning: Only triclinic boxes with the first vector parallel to the x-axis
and the second vector in the xy-plane are supported.
Box (3x3):
Box[0]={ nan, 0.0e+00, 0.0e+00}
Box[1]={ nan, nan, nan}
Box[2]={ nan, nan, nan}
Can not fix pbc.
;Run control
integrator = md
dt = 0.002
nsteps = 100 ;5000
nstcomm = 100
;Output control
nstenergy = 100
nstxout = 100
nstvout = 0
nstfout = 0
nstlog = 1000
nstxtcout = 1000
;Neighbor searching
nstlist = 10
ns_type = grid
;Electrostatics/VdW
coulombtype = Shift
vdw-type= Shift
rcoulomb-switch = 0
rvdw-switch = 0.9 ;0
;Cut-offs
rlist = 1.25
rcoulomb= 1.0
rvdw= 1.0
;Temperature coupling
Tcoupl = v-rescale
tc-grps = System
tau_t = 0.1
ref_t = 300
;Pressure coupling
Pcoupl = Parrinello-Rahman
Pcoupltype = isotropic
tau_p = 1
compressibility = 3.5e-5
ref_p = 10
;Velocity generation
gen_vel = yes
gen_temp= 300.0
gen_seed= 173529
;Bonds
constraints = all-bonds
constraint-algorithm = lincs
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www interface
or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] error: Only triclinic boxes...
Hello all,
I am getting the error below at the very beginning of the simulation (both
serial and parallel). I am sure I did not encounter this problem before with
the same input files. This has just happened now. I really have no clue why
this is happening. could you please help me? Thank you all in advance.
Warning: Only triclinic boxes with the first vector parallel to the x-axis
and the second vector in the xy-plane are supported.
Box (3x3):
Box[0]={ nan, 0.0e+00, 0.0e+00}
Box[1]={ nan, nan, nan}
Box[2]={ nan, nan, nan}
Can not fix pbc.
;Run control
integrator = md
dt = 0.002
nsteps = 100 ;5000
nstcomm = 100
;Output control
nstenergy = 100
nstxout = 100
nstvout = 0
nstfout = 0
nstlog = 1000
nstxtcout = 1000
;Neighbor searching
nstlist = 10
ns_type = grid
;Electrostatics/VdW
coulombtype = Shift
vdw-type= Shift
rcoulomb-switch = 0
rvdw-switch = 0.9 ;0
;Cut-offs
rlist = 1.25
rcoulomb= 1.0
rvdw= 1.0
;Temperature coupling
Tcoupl = v-rescale
tc-grps = System
tau_t = 0.1
ref_t = 300
;Pressure coupling
Pcoupl = Parrinello-Rahman
Pcoupltype = isotropic
tau_p = 1
compressibility = 3.5e-5
ref_p = 10
;Velocity generation
gen_vel = yes
gen_temp= 300.0
gen_seed= 173529
;Bonds
constraints = all-bonds
constraint-algorithm = lincs
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] error: Only triclinic boxes...
Hello Justin,
Thank you. I am using
gmxdump -s *.tpr -om to produce mdout.mdp file but I dont see box vector
information. This error happens at the very beginning. I see no steps are
calculated.
On 2 June 2011 13:11, Justin A. Lemkul jalem...@vt.edu wrote:
Elisabeth wrote:
Hello all,
I am getting the error below at the very beginning of the simulation (both
serial and parallel). I am sure I did not encounter this problem before with
the same input files. This has just happened now. I really have no clue why
this is happening. could you please help me? Thank you all in advance.
Warning: Only triclinic boxes with the first vector parallel to the x-axis
and the second vector in the xy-plane are supported.
Box (3x3):
Box[0]={ nan, 0.0e+00, 0.0e+00}
Box[1]={ nan, nan, nan}
Box[2]={ nan, nan, nan}
Can not fix pbc.
Your system is probably blowing up. Nan means not a number, so the box
vectors have become either infinitely large or small. Either the input
coordinate file contained a malformed or non-existent box, or the simulation
is collapsing along the way somewhere. Does the simulation run for a while
before printing this, or is it right away?
You can check the starting box vectors in the .tpr file with gmxdump to
verify that they are sensible.
-Justin
;Run control integrator = md
dt = 0.002 nsteps =
100 ;5000 nstcomm = 100
;Output control
nstenergy = 100nstxout = 100
nstvout = 0
nstfout = 0
nstlog = 1000 nstxtcout = 1000
;Neighbor searching
nstlist = 10 ns_type = grid
;Electrostatics/VdW
coulombtype = Shiftvdw-type=
Shift rcoulomb-switch = 0 rvdw-switch
= 0.9 ;0
;Cut-offs
rlist = 1.25 rcoulomb= 1.0
rvdw= 1.0
;Temperature coupling Tcoupl = v-rescale
tc-grps = System tau_t = 0.1
ref_t = 300 ;Pressure coupling
Pcoupl = Parrinello-RahmanPcoupltype
= isotropic tau_p = 1
compressibility = 3.5e-5 ref_p = 10
;Velocity generation gen_vel = yes
gen_temp= 300.0 gen_seed=
173529
;Bonds
constraints = all-bonds constraint-algorithm =
lincs
--
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www interface
or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] error: Only triclinic boxes...
Dear Justin,
I find:
box (3x3):
box[0]={ 1.6e+01, 0.0e+00, 0.0e+00}
box[1]={ 0.0e+00, 1.6e+01, 0.0e+00}
box[2]={ 0.0e+00, 0.0e+00, 1.6e+01}
box_rel (3x3):
box_rel[0]={ 0.0e+00, 0.0e+00, 0.0e+00}
box_rel[1]={ 0.0e+00, 1.0e+00, 0.0e+00}
box_rel[2]={ 0.0e+00, 0.0e+00, 1.0e+00}
boxv (3x3):
boxv[0]={ 0.0e+00, 0.0e+00, 0.0e+00}
boxv[1]={ 0.0e+00, 0.0e+00, 0.0e+00}
boxv[2]={ 0.0e+00, 0.0e+00, 0.0e+00}
pres_prev (3x3):
pres_prev[0]={ 0.0e+00, 0.0e+00, 0.0e+00}
pres_prev[1]={ 0.0e+00, 0.0e+00, 0.0e+00}
pres_prev[2]={ 0.0e+00, 0.0e+00, 0.0e+00}
svir_prev (3x3):
svir_prev[0]={ 0.0e+00, 0.0e+00, 0.0e+00}
svir_prev[1]={ 0.0e+00, 0.0e+00, 0.0e+00}
svir_prev[2]={ 0.0e+00, 0.0e+00, 0.0e+00}
fvir_prev (3x3):
fvir_prev[0]={ 0.0e+00, 0.0e+00, 0.0e+00}
fvir_prev[1]={ 0.0e+00, 0.0e+00, 0.0e+00}
fvir_prev[2]={ 0.0e+00, 0.0e+00, 0.0e+00}
nosehoover_xi: not available
x (2896x3):
x[0]={-7.49400e+00, 4.57000e-01, 5.0e-01}
x[1]={-7.58200e+00, 5.21000e-01, 5.0e-01}
x[2]={-7.49600e+00, 3.94000e-01, 4.11000e-01}
x[3]={-7.49600e+00, 3.94000e-01, 5.89000e-01}
x[4]={-7.36800e+00, 5.43000e-01, 5.0e-01}
x[5]={-7.36800e+00, 6.06000e-01, 5.89000e-01}
x[6]={-7.36800e+00, 6.06000e-01, 4.11000e-01}
x[7]={-7.24200e+00, 4.57000e-01, 5.0e-01}
x[8]={-7.24200e+00, 3.94000e-01, 4.11000e-01}
x[9]={-7.24200e+00, 3.94000e-01, 5.89000e-01}
x[ 10]={-7.11700e+00, 5.43000e-01, 5.0e-01}
x[ 11]={-7.11700e+00, 6.06000e-01, 5.89000e-01}
x[ 12]={-7.11700e+00, 6.06000e-01, 4.11000e-01}
x[ 13]={-6.99100e+00, 4.57000e-01, 5.0e-01}
I made box bigger to remove bad contacts in initial configuration (although
I have been using the same gro file before and never got such error). Do you
recommend I reinstall 4.5.4 ?
Appreciate any help...
Regards,
On 2 June 2011 14:47, Justin A. Lemkul jalem...@vt.edu wrote:
Elisabeth wrote:
Hello Justin,
Thank you. I am using
gmxdump -s *.tpr -om to produce mdout.mdp file but I dont see box
You don't need an mdout.mdp file, nor would box vectors be present in one.
What you're looking for is the box information present in the .tpr file.
Run:
gmxdump -s topol.tpr out
Then search in the out file for the box information.
vector information. This error happens at the very beginning. I see no
steps are calculated.
Then the input box information is likely corrupt, or the simulation blows
up at the very start of the simulation (i.e. step 0).
-Justin
On 2 June 2011 13:11, Justin A. Lemkul jalem...@vt.edu mailto:
jalem...@vt.edu wrote:
Elisabeth wrote:
Hello all,
I am getting the error below at the very beginning of the
simulation (both serial and parallel). I am sure I did not
encounter this problem before with the same input files. This
has just happened now. I really have no clue why this is
happening. could you please help me? Thank you all in advance.
Warning: Only triclinic boxes with the first vector parallel to
the x-axis and the second vector in the xy-plane are supported.
Box (3x3):
Box[0]={ nan, 0.0e+00, 0.0e+00}
Box[1]={ nan, nan, nan}
Box[2]={ nan, nan, nan}
Can not fix pbc.
Your system is probably blowing up. Nan means not a number, so
the box vectors have become either infinitely large or small.
Either the input coordinate file contained a malformed or
non-existent box, or the simulation is collapsing along the way
somewhere. Does the simulation run for a while before printing
this, or is it right away?
You can check the starting box vectors in the .tpr file with gmxdump
to verify that they are sensible.
-Justin
;Run control integrator
= mddt = 0.002
nsteps = 100 ;5000 nstcomm=
100 ;Output control
nstenergy = 100nstxout =
100 nstvout = 0
nstfout = 0
nstlog = 1000 nstxtcout = 1000
;Neighbor searching
nstlist = 10 ns_type =
grid ;Electrostatics/VdW
coulombtype
[gmx-users] compression to study pressure effect
Dear all, I am trying to study the effect of pressure on total potential of my system (8 polymer chains). My problem is that I dont see a systematic effect of pressure on potentials and I cant judge if different pressures increase or decrease potential. This is the critical observable in my system and with high fluctuations I am getting I cant make comment on pressure effect. Total drift is high in most potential functions..I start my runs from a frame which is the output of an older NPT run (I use cpt file) that has a close density to what I want. In the production runs (NPT for 12 27 70 bar). Below is the settings I am using and I really appreciate it if you could comment on the most important factors for such a study. I tried different cutoffs as well...I though maybe increasing cutoffs invloves more interactions and can represent better the pressure effect... Thanks so much! ;Run control integrator = md dt = 0.002 nsteps = 100 ;5000 nstcomm = 100 ;Output control nstenergy = 100 nstxout = 100 nstvout = 0 nstfout = 0 nstlog = 1000 nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0.5 ;0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist = 1.32 ;1.25 ; tired different cutoffs and r -switch... rcoulomb= 1.1 ;1.0 rvdw = 1.1 ;1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 ;Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 compressibility = 3.5e-5 ref_p = 10 ;Velocity generation gen_vel = no;yes gen_temp= 300.0 gen_seed= 173529 ;Bonds constraints = all-bonds constraint-algorithm = lincs results for 12 27 70 bar resp. Statistics over 2043101 steps [ 500. through 4586.2000 ps ], 17 data sets All statistics are over 204311 points Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-634.3582.911.0696 -20.3551 (kJ/mol) Coulomb (SR)184.0821.43.627129.04046 (kJ/mol) P*otential 642.5286.423.1897 -44.6636 (kJ/mol)* Kinetic En. 904.294 0.07516.2742 0.382223 (kJ/mol) Total Energy1546.826.328.6865 -44.2813 (kJ/mol) Temperature 300.186 0.025 5.4023 0.126881 (K) Pressure12.3856 0.89887.5015.41599 (bar) Density 824.4172.19.2783814.7965 (kg/m^3) Statistics over 750001 steps [ 500. through 2000. ps ], 17 data sets All statistics are over 75001 points Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-638.6434.312.0718 -10.2094 (kJ/mol) Coulomb (SR) 184.51 12.897775.99117 (kJ/mol) *Potential 634.9025.722.2857 -21.4761 (kJ/mol)* Kinetic En. 904.181 0.1116.2261 -0.275076 (kJ/mol) Total Energy1539.085.727.9661 -21.7511 (kJ/mol) Temperature 300.148 0.0355.38636 -0.0913135 (K) Pressure27.64671.4881.4098.58275 (bar) Density 826.7343.49.760896.06824 (kg/m^3) 70 bar run even 6 ns runs gives large total drifts..e.g.13.4934 for density.. Statistics over 2750001 steps [ 500. through 6000. ps ], 17 data sets All statistics are over 275001 points Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-639.2063.711.6741 -23.7773 (kJ/mol) Coulomb (SR) 182.971.74.1866311.8694 (kJ/mol) *Potential 640.1147.524.9191 -51.8187 (kJ/mol)* Kinetic En. 904.3 0.07416.2887 0.161695 (kJ/mol) Total Energy1544.417.530.1347-51.657 (kJ/mol) Temperature 300.187 0.0245.40713 0.0536754 (K) Pressure69.7658 0.86893.933 -2.29885 (bar) Density 829.6572.38.9884913.4934 (kg/m^3) --
[gmx-users] compressing with NPT
Dear experts, I have been considering effect of fixing bond lengths in my system (constraints= none / all-bonds) on equilibrated density using NPT runs with 4.5.4. I am showing two set of results for two ref pressures of 10 and 70 bar. in case of constraints= none I had to employ a 0.5 fs timestep otherwise system blows up. Can you please comment why equilibrated densities are different for the same reference pressure? Thanks. 70 bar, all bonds constraint timestep 2fs Statistics over 875001 steps [ 250. through 2000. ps ], 2 data sets All statistics are over 87501 points Energy Average Err.Est. RMSD Tot-Drift --- Pressure70.4993 0.092542.944 0.222707 (bar) Density 625.894 0.3911.9075 0.110708 (kg/m^3) 70 bar, constraints= none timestep 0.5 fs Statistics over 351 steps [ 250. through 2000. ps ], 2 data sets All statistics are over 350001 points Energy Average Err.Est. RMSD Tot-Drift --- Pressure70.6207 0.11932.219 -0.500892 (bar) Density 614.451 0.5513.7381 -0.347349 (kg/m^3) --- 10 bar, all bonds constraint timestep 2fs Statistics over 875001 steps [ 250. through 2000. ps ], 2 data sets All statistics are over 87501 points Energy Average Err.Est. RMSD Tot-Drift --- Pressure10.4843 0.039621.684 0.0156492 (bar) Density 624.362 112.1354 -2.11725 (kg/m^3) 10 bar, constraints= none timestep 0.5 fs Statistics over 51 steps [ 250. through 500. ps ], 2 data sets All statistics are over 50001 points Energy Average Err.Est. RMSD Tot-Drift --- Temperature 299.92 0.114.86352 -0.139097 (K) Density 600.689 214.759711.6485 (kg/m^3) ;Run control integrator = md dt = 0.002 ; 0.0005 for constraints =none otherwise blows up nsteps = 100 nstcomm = 100 ;Output control nstenergy = 100 nstxout = 100 nstvout = 0 nstfout = 0 nstlog = 1000 nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 ;Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 2 compressibility = 4.5e-5 4.5e-5 ref_p = 10 constraints = all-bonds constraint-algorithm = lincs -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] NPT for compressing
On 20 May 2011 23:04, Mark Abraham mark.abra...@anu.edu.au wrote: On 21/05/2011 10:45 AM, Elisabeth wrote: Dear experts, I have a box of pure hexane with density of around 50 SI (size 7nm). With the settings below I get error: vol 0.65 imb F 3% step 42100, will finish Fri May 20 18:54:05 2011 Warning: 1-4 interaction between 1246 and 1249 at distance 3.580 which is larger than the 1-4 table size 2.250 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Warning: 1-4 interaction between 1249 and 1252 at distance 3.693 which is larger than the 1-4 table size 2.250 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size A list of missing interactions: Bond of 2375 missing 2 Angle of 4500 missing 11 Ryckaert-Bell. of 5625 missing 21 LJ-14 of 5625 missing 15 Molecule type 'Hexane' the first 10 missing interactions, except for exclusions: Ryckaert-Bell. atoms58 11 14 global 1245 1248 1251 1254 LJ-14 atoms5 14 global 1245 1254 Angle atoms8 11 14 global 1248 1251 1254 Ryckaert-Bell. atoms8 11 14 15 global 1248 1251 1254 1255 Ryckaert-Bell. atoms8 11 14 16 global 1248 1251 1254 1256 Ryckaert-Bell. atoms8 11 14 17 global 1248 1251 1254 1257 LJ-14 atoms8 15 global 1248 1255 LJ-14 atoms8 16 global 1248 1256 Ryckaert-Bell. atoms98 11 14 global 1249 1248 1251 1254 LJ-14 atoms9 14 global 1249 1254 --- Program mdrun_mpi, VERSION 4.5.4 Source code file: domdec_top.c, line: 356 Fatal error: 49 of the 18125 bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (1.25 nm) or the two-body cut-off distance (1.25 nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors OK, your system blew up. mdp file: ;Bonds constraints = none constraint-algorithm = lincs ;Run control integrator = md dt = 0.001 nsteps = 50 ;5000 nstcomm = 100 ;Output control nstenergy = 100 ; frequency to write energies to energy file. i.e., energies and other statistical data are stored every 10 steps nstxout = 100;1 nstvout = 0 nstfout = 0 nstlog = 1000 ; frequency to write energies to log file nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 ;Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 compressibility = 4.5e-5 4.5e-5 ref_p = 10 ;Velocity generation gen_vel = yes gen_temp= 300.0 gen_seed= 173529 Generating velocities and compressing from an unequilibrated initial conformation in the same run is asking for trouble, and sometimes you're observing it. See http://www.gromacs.org/Documentation/How-tos/Steps_to_Perform_a_Simulation. Equilibrate first. Then slowly change the conditions so that you're never all that far from equilibrium, and reach your destination. Hello Mark, My box of alkane was already equilibrated using a NVT run. I just have to increase the density. After NVT run, do I have to turn off gen_vel = yes for NPT run? Is this necessary for all NPT runs? When I use constraints for all-bonds and dt of 0.002 compressing the box works well and I get density of above 600 SI which is the actual density of alkane. Can any one help why setting no constraints and dt=0.001 leads to the error message shown above. You probably got lucky the first time and not the second. Or, your atom-atom bond interactions are no good. When
Re: [gmx-users] NPT for compressing
On 21 May 2011 19:02, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: On 20 May 2011 23:04, Mark Abraham mark.abra...@anu.edu.au mailto: mark.abra...@anu.edu.au wrote: On 21/05/2011 10:45 AM, Elisabeth wrote: Dear experts, I have a box of pure hexane with density of around 50 SI (size 7nm). With the settings below I get error: vol 0.65 imb F 3% step 42100, will finish Fri May 20 18:54:05 2011 Warning: 1-4 interaction between 1246 and 1249 at distance 3.580 which is larger than the 1-4 table size 2.250 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Warning: 1-4 interaction between 1249 and 1252 at distance 3.693 which is larger than the 1-4 table size 2.250 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size A list of missing interactions: Bond of 2375 missing 2 Angle of 4500 missing 11 Ryckaert-Bell. of 5625 missing 21 LJ-14 of 5625 missing 15 Molecule type 'Hexane' the first 10 missing interactions, except for exclusions: Ryckaert-Bell. atoms58 11 14 global 1245 1248 1251 1254 LJ-14 atoms5 14 global 1245 1254 Angle atoms8 11 14 global 1248 1251 1254 Ryckaert-Bell. atoms8 11 14 15 global 1248 1251 1254 1255 Ryckaert-Bell. atoms8 11 14 16 global 1248 1251 1254 1256 Ryckaert-Bell. atoms8 11 14 17 global 1248 1251 1254 1257 LJ-14 atoms8 15 global 1248 1255 LJ-14 atoms8 16 global 1248 1256 Ryckaert-Bell. atoms98 11 14 global 1249 1248 1251 1254 LJ-14 atoms9 14 global 1249 1254 --- Program mdrun_mpi, VERSION 4.5.4 Source code file: domdec_top.c, line: 356 Fatal error: 49 of the 18125 bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (1.25 nm) or the two-body cut-off distance (1.25 nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors OK, your system blew up. mdp file: ;Bonds constraints = noneconstraint-algorithm = lincs ;Run control integrator = md dt = 0.001 nsteps = 50 ;5000 nstcomm = 100 ;Output control nstenergy = 100 ; frequency to write energies to energy file. i.e., energies and other statistical data are stored every 10 steps nstxout = 100;1nstvout = 0 nstfout = 0 nstlog = 1000 ; frequency to write energies to log file nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = Shift vdw-type = Shiftrcoulomb-switch = 0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist = 1.25 rcoulomb= 1.0rvdw= 1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 ;Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 compressibility = 4.5e-5 4.5e-5 ref_p = 10 ;Velocity generation gen_vel = yesgen_temp= 300.0 gen_seed= 173529 Generating velocities and compressing from an unequilibrated initial conformation in the same run is asking for trouble, and sometimes you're observing it. See http://www.gromacs.org/Documentation/How-tos/Steps_to_Perform_a_Simulation . Equilibrate first. Then slowly change the conditions so that you're never all that far from equilibrium, and reach your
[gmx-users] NPT for compressing
Dear experts, I have a box of pure hexane with density of around 50 SI (size 7nm). With the settings below I get error: vol 0.65 imb F 3% step 42100, will finish Fri May 20 18:54:05 2011 Warning: 1-4 interaction between 1246 and 1249 at distance 3.580 which is larger than the 1-4 table size 2.250 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Warning: 1-4 interaction between 1249 and 1252 at distance 3.693 which is larger than the 1-4 table size 2.250 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size A list of missing interactions: Bond of 2375 missing 2 Angle of 4500 missing 11 Ryckaert-Bell. of 5625 missing 21 LJ-14 of 5625 missing 15 Molecule type 'Hexane' the first 10 missing interactions, except for exclusions: Ryckaert-Bell. atoms58 11 14 global 1245 1248 1251 1254 LJ-14 atoms5 14 global 1245 1254 Angle atoms8 11 14 global 1248 1251 1254 Ryckaert-Bell. atoms8 11 14 15 global 1248 1251 1254 1255 Ryckaert-Bell. atoms8 11 14 16 global 1248 1251 1254 1256 Ryckaert-Bell. atoms8 11 14 17 global 1248 1251 1254 1257 LJ-14 atoms8 15 global 1248 1255 LJ-14 atoms8 16 global 1248 1256 Ryckaert-Bell. atoms98 11 14 global 1249 1248 1251 1254 LJ-14 atoms9 14 global 1249 1254 --- Program mdrun_mpi, VERSION 4.5.4 Source code file: domdec_top.c, line: 356 Fatal error: 49 of the 18125 bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (1.25 nm) or the two-body cut-off distance (1.25 nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors mdp file: ;Bonds constraints = none constraint-algorithm = lincs ;Run control integrator = md dt = 0.001 nsteps = 50 ;5000 nstcomm = 100 ;Output control nstenergy = 100 ; frequency to write energies to energy file. i.e., energies and other statistical data are stored every 10 steps nstxout = 100;1 nstvout = 0 nstfout = 0 nstlog = 1000 ; frequency to write energies to log file nstxtcout = 1000 ;Neighbor searching nstlist = 10 ns_type = grid ;Electrostatics/VdW coulombtype = Shift vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 ;0 ;Cut-offs rlist = 1.25 rcoulomb= 1.0 rvdw= 1.0 ;Temperature coupling Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 ;Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 compressibility = 4.5e-5 4.5e-5 ref_p = 10 ;Velocity generation gen_vel = yes gen_temp= 300.0 gen_seed= 173529 When I use constraints for all-bonds and dt of 0.002 compressing the box works well and I get density of above 600 SI which is the actual density of alkane. Can any one help why setting no constraints and dt=0.001 leads to the error message shown above. When no constraints is used PR doesnt work for compressing and berendsen gives desired densty only for high pressures above 50 bar. I need to study my system at different pressures from 10 bar... please comment on this.. Thanks. 10 bar ; Pressure coupling Pcoupl = berendsen Pcoupltype = isotropic tau_p = * 0.1 * 0.5 compressibility = 4.5e-5 4.5e-5 ref_p = 10 10 does not compress! With and without constr. Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] How to freeze groups
Hello, Can anyone suggest some tutorials or sample top files for simulations involving freezing of groups? I dont see much info in chapter 5 of manual. Thanks, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] multicomponent system- units
Hi Mark, I am excited to see that there is a solution to my issue. I thought this problem can not be resolved. In thermodynamics of polymer solutions, people use some models (equation of state) in which an interaction parameter K_AB appears which is defined in terms of interaction energies i.e. 1-K_AB=(E_AB)/(E_AA*E_BB)^0.5. One way to obtain this parameter is to manipulate this K so that equation of state predicts say bubble point data or density vs. pressure. In this procedure they dont look at interaction energies E_BB,...and only K is tuned. (or in some models they deal with E_ij interaction energies and manipulate so that some properties are fitted to experimental data). Now what I am interested in is calculating these interaction energies by MD and thats why I need to extract pairwise energies per mol. To double check what I have done with you: FOr a system having 4 polymer chains and 100 solvent molecules, I defined two groups in index file: [polymer] with all atoms of polymer chains. and [solvent] with all atoms of solvent. and use energygrps= polymer solvent. Now I have polymer-solvent, polymer-polymer and solvent-solvent interaction energies (LJ + Coulomb SR for each pair). As you say to normalize this I have to divide by [(4*Np)*(100*Ns)] where Np and Ns are number of atoms in polymer chain and solvent molecule. 1- Did I get your instruction correctly? 2- The unit of energies is per atom now? I am confused if its per atom or molecule? 3- Since the interaction parameter in the model is defined as 1- K_AB=(E_AB)/(E_AA*E_BB)^0.5 and the ratio of interaction energies appear in K, is this normalization sufficient? I mean because of ratio of energies it seems there is no need to convert these normalized values to MOL! 4- Is it possible to achieve energy per MOL for this binary system from normalized energies? Appreciate your help! Best :) On 12 April 2011 00:10, Mark Abraham mark.abra...@anu.edu.au wrote: Hello Mark, Thank you for your reply. I have already created the energy groups. I am trying to validate pairwise energy values (nonbonded) with some other work ( a thermodynamic model) where they fit these AA AB BB (E_AA, E_AB, E_BB) energies so that some phase diagrams are reproduced. The pairwise energies defined in the model are in KJ/mol. So how did they compute these interaction energies? The energy quantity GROMACS reports for a microstate can be best thought of as the energy one would have for a mole of such microstates. Alternatively, divide by N_A and that's the energy for this microstate - but that's a much less convenient number to use. To obtain a quantity that is independent of the number of particles, you have to normalize for the number of interactions of each type. If these are all pairwise between atoms in a unary system, then you need to divide by the square of the number of atoms. So for the mixed interaction energy of the binary system, you divide by the product of the respective numbers of atoms. You should also verify that these actually are converged observables that are independent of the number of particles by simulating replicates from different starting configurations, and systems of different sizes. Mark Since my energies are not per mol, my results are useless, unfortunately. As they depend on number of molecules in the system. To achieve my goal, what do you suggest? For a binary system, can I run two separate simulations for pure A and B in which case using -nmol gives per mol energies and somehow predict AB from them? Does this make sense? Please guide me, I am stuck on this.. Thanks, On 9 April 2011 20:56, Mark Abraham mark.abra...@anu.edu.au wrote: On 8/04/2011 12:18 PM, Elisabeth wrote: Hello everyone, I have encountered a simple problem. For a homogenous system what g_energy reports is dependent on the system size and one needs to use -nmol option to divide energies by number of molecules to obtain per mol values. I am attempting to extract interaction energies between species in a three component system. I am puzzled how this can be achieved for such a system. Say there are 100 solvent, 20 solute A and 10 B molecules. You would have to start by defining energy groups that contain relevant sets of molecules (see manual). Even once you've got them, the group-wise energies won't mean much of anything. Every observable is dependent on the configuration microstate, and unless you can estimate the relative population of different microstates to estimate a free energy... Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing
Re: [gmx-users] multicomponent system- units
Hello Mark, Thank you for your reply. I have already created the energy groups. I am trying to validate pairwise energy values (nonbonded) with some other work ( a thermodynamic model) where they fit these AA AB BB (E_AA, E_AB, E_BB) energies so that some phase diagrams are reproduced. The pairwise energies defined in the model are in KJ/mol. Since my energies are not per mol, my results are useless, unfortunately. As they depend on number of molecules in the system. To achieve my goal, what do you suggest? For a binary system, can I run two separate simulations for pure A and B in which case using -nmol gives per mol energies and somehow predict AB from them? Does this make sense? Please guide me, I am stuck on this.. Thanks, On 9 April 2011 20:56, Mark Abraham mark.abra...@anu.edu.au wrote: On 8/04/2011 12:18 PM, Elisabeth wrote: Hello everyone, I have encountered a simple problem. For a homogenous system what g_energy reports is dependent on the system size and one needs to use -nmol option to divide energies by number of molecules to obtain per mol values. I am attempting to extract interaction energies between species in a three component system. I am puzzled how this can be achieved for such a system. Say there are 100 solvent, 20 solute A and 10 B molecules. You would have to start by defining energy groups that contain relevant sets of molecules (see manual). Even once you've got them, the group-wise energies won't mean much of anything. Every observable is dependent on the configuration microstate, and unless you can estimate the relative population of different microstates to estimate a free energy... Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] multicomponent system- units
Hello everyone, I have encountered a simple problem. For a homogenous system what g_energy reports is dependent on the system size and one needs to use -nmol option to divide energies by number of molecules to obtain per mol values. I am attempting to extract interaction energies between species in a three component system. I am puzzled how this can be achieved for such a system. Say there are 100 solvent, 20 solute A and 10 B molecules. Thanks! Regards -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] PME
Hello Justin, Several days ago you answered my question about calculating nonbonded terms: Question: If I want to look at nonboded interactions only, do I have to add Coul. recip. to [ LJ (SR) + Coulomb (SR) ] ? Answer: The PME-related terms contain both solute-solvent, solvent-solvent, and potentially solute-solute terms (depending on the size and nature of the solute), so trying to interpret this term in some pairwise fashion is an exercise in futility. my question is if I want to add up nonbonded related terms to get inter molecular energies, do I have to add Coul. recip. or it is already included in Coulomb (SR)? and also, for a A-B system, I have been using energy groups to extract solute-solvent, solvent-solvent, solute-solute terms. Did you mean that applying doing so with PME as electrostatics treatment is not correct? Thanks for your help! Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] PME
On 6 April 2011 15:01, Michael Brunsteiner mbx0...@yahoo.com wrote: Elisabeth, You CAN, in fact calculate the contribution of the reciprocal part of the PME energy to the binding energy between two components in a heterogeneous system, its just quite tedious... say, your system is molecules A and B for which you want to know the interaction energy, and the rest of the system, typically the solvent, we call C. Now your total Reciprocal Coulomb energy will have six parts: ER_tot = ER_AA + ER_BB + ER_CC + ER_AB + ER_AC + ER_BC but these parts are NOT given in the gromacs output as they cannot be calculated DIRECTLY, you have to calculate them by setting the charges on A, B, or C (or combinations thereof) to zero (there is a tool for setting the charges in a tpr file to zero) and then do more runs with: mdrun -rerun based on the original trajectory to get the required contributions. then E_AB = ER_C0 - ER_A0C0 - ER_B0C0 (or something like it, do double check that formula, i can't be bothered thinking it through now ... here ER_A0C0, for example, is the reciprocal part of the coulomb energy with charges in groups A and C set to zero, etc) this being said ... it's tedious, time-consuming, and error-prone (you need to use double precision and save a lot of frames to get reasonably accurate numbers) You might be better off using reaction field, or PME and simply ignore the reciprocal part altogether (if your molecules A, B are NOT charged and have no permanent and large dipole moment you might get away with the latter) Thanks for your elaborate message. The point is in my case there is no option other than ignoring LR since LR is not covered by shift or switch functions but at least what PME reports for SR is more accurate. So the decomposed Coulmb. SR terms I am getting using energy groups from PME are reliable ? BTW: I am dealing with non polar particles i.e alkanes and carbon and hydrogen are the only species I have. Can you please tell me about the tool in tpr file that sets all charges to zero..I might use this to check how turning off electrostatics affects properties. and just a little question: I am unclear about LJ-14 and Coulomb-14 too. Are these included in LJ-SR and Coulomb-SR or for each pair one needs to add up the respective 14 term? i.e A-B LJ-14 + A-B LJ-SR + A-B Coulomb-14 + A-B Coulomb-SR to get nonbonded inter molecular energy for A-B components? If they are already included what is the point of reporting them separately? Thank you so much, What Justin said is correct, PME (or any other Ewald-like method, PPPM, FMA, etc) is standard these days, and for a good reason. However, different properties are affected to a different extent by neglecting the long range interactions, and for what you want to calculate it might be OK for getting at least a qualitative answer, as long as you use PME for the actual MD. (I'd be VERY surprised if everybody who did LIE in the last 10 years went through the trouble outlined above) have fun! mic Elisabeth wrote: Hello Justin, Several days ago you answered my question about calculating nonbonded terms: Question: If I want to look at nonboded interactions only, do I have to add Coul. recip. to [ LJ (SR) + Coulomb (SR) ] ? Answer: The PME-related terms contain both solute-solvent, solvent-solvent, and potentially solute-solute terms (depending on the size and nature of the solute), so trying to interpret this term in some pairwise fashion is an exercise in futility. my question is if I want to add up nonbonded related terms to get inter molecular energies, do I have to add Coul. recip. or it is already included in Coulomb (SR)? They are separate energy terms. The PME mesh terms is Coul. recip. and the short-range interactions (contained within rcoulomb, calculated by a modified switch potential) are Coulomb (SR). and also, for a A-B system, I have been using energy groups to extract solute-solvent, solvent-solvent, solute-solute terms. Did you mean that applying doing so with PME as electrostatics treatment is not correct? PME has been consistently shown to be one of the most accurate long-range electrostatics methods and is widely used, but in your case is preventing you from extracting the quantity you're after (if it can even be reasonably defined at all). Using energygrps will not resolve the problem I described above. The Coul. recip. term contains long-range energies between (potentially) A-B, A-A, and B-B, depending on the nature of what A and B are. The only terms that are decomposed via energygrps are the short-range terms, which are calculated pairwise. Thus, with PME, there is no straightforward way to simply define an intermolecular energy for a heterogeneous system. You might be able to define such a term for a completely homogeneous system (which also assumes that the sampling has converged
Re: [gmx-users] PME
On 6 April 2011 19:28, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: On 6 April 2011 15:01, Michael Brunsteiner mbx0...@yahoo.com mailto: mbx0...@yahoo.com wrote: Elisabeth, You CAN, in fact calculate the contribution of the reciprocal part of the PME energy to the binding energy between two components in a heterogeneous system, its just quite tedious... say, your system is molecules A and B for which you want to know the interaction energy, and the rest of the system, typically the solvent, we call C. Now your total Reciprocal Coulomb energy will have six parts: ER_tot = ER_AA + ER_BB + ER_CC + ER_AB + ER_AC + ER_BC but these parts are NOT given in the gromacs output as they cannot be calculated DIRECTLY, you have to calculate them by setting the charges on A, B, or C (or combinations thereof) to zero (there is a tool for setting the charges in a tpr file to zero) and then do more runs with: mdrun -rerun based on the original trajectory to get the required contributions. then E_AB = ER_C0 - ER_A0C0 - ER_B0C0 (or something like it, do double check that formula, i can't be bothered thinking it through now ... here ER_A0C0, for example, is the reciprocal part of the coulomb energy with charges in groups A and C set to zero, etc) this being said ... it's tedious, time-consuming, and error-prone (you need to use double precision and save a lot of frames to get reasonably accurate numbers) You might be better off using reaction field, or PME and simply ignore the reciprocal part altogether (if your molecules A, B are NOT charged and have no permanent and large dipole moment you might get away with the latter) Thanks for your elaborate message. The point is in my case there is no option other than ignoring LR since LR is not covered by shift or switch functions but at least what PME reports for SR is more accurate. So the decomposed Coulmb. SR terms I am getting using energy groups from PME are reliable ? I don't understand your question entirely, so hopefully someone else can comment. Hi Justin, I am using PME and extract decomposed Coulmb. SR terms using energy groups from g_energy. As we discussed LR terms (coulmb recip) can not be decomposed. What I want to make sure about is that at least energy groups give reliable PME Coulmb. SR terms.. Reading your statement below makes me interpret that both PME related terms i.e SR and LR (coulmb recip.) can no be decomposed. so again I am copying your statement : The *PME-related terms* contain both solute-solvent, solvent-solvent, and potentially solute-solute terms (depending on the size and nature of the solute), so trying to interpret this term in some pairwise fashion is an exercise in futility. In other words *if one needs to obtain decomposed nonbonded intermolecular terms*, PME is not an option and maybe shift potentials must be used. Is that what you mean? I appreciate any clarification on *PME-related terms*... Thanks :) Best, BTW: I am dealing with non polar particles i.e alkanes and carbon and hydrogen are the only species I have. Can you please tell me about the tool in tpr file that sets all charges to zero..I might use this to check how turning off electrostatics affects properties. tpbconv -zeroq and just a little question: I am unclear about LJ-14 and Coulomb-14 too. Are these included in LJ-SR and Coulomb-SR or for each pair one needs to add up the respective 14 term? i.e A-B LJ-14 + A-B LJ-SR + A-B Coulomb-14 + A-B Coulomb-SR to get nonbonded inter molecular energy for A-B components? If they are already included what is the point of reporting them separately? 1-4 interactions are intramolecular, not intermolecular. Every nonbonded energy term that is listed in the .edr file is a separate entity. -Justin Thank you so much, What Justin said is correct, PME (or any other Ewald-like method, PPPM, FMA, etc) is standard these days, and for a good reason. However, different properties are affected to a different extent by neglecting the long range interactions, and for what you want to calculate it might be OK for getting at least a qualitative answer, as long as you use PME for the actual MD. (I'd be VERY surprised if everybody who did LIE in the last 10 years went through the trouble outlined above) have fun! mic Elisabeth wrote: Hello Justin, Several days ago you answered my question about calculating nonbonded terms: Question: If I want to look at nonboded interactions only, do I have to add Coul. recip. to [ LJ (SR) + Coulomb (SR) ] ? Answer: The PME-related terms contain both solute-solvent, solvent-solvent, and potentially solute-solute terms (depending on the size and nature of the solute), so trying to interpret
Re: [gmx-users] Heat of vap
Elisabeth wrote: Dear David, I followed your instructions and calculated Heat of vaporization of my alkane once with one molecule in gas phase (no cutoff) and once with equivalent number of molecules as in liquid phase as Justin suggested. Results are as follows: To get heat of vaporization, you shouldn't be simulating just a single molecule in the gas phase, it should be an equivalent number of molecules as you have in the liquid phase. Hello David and Justin, My explanation was not clear. Below is the results for liquid phase and for gas phase I tried two cases: one single molecule and the other time for equivalent number of molecules as in liquid phase and thats why results are very similar. ( However Justin says one single molecule is not correct. I think when cutoffs is set to zero only bonded terms are treated and even where there are many particles in gas phase to get energies per mole of molecules i.e g_energy -nmol XXX must be used so values should be colse to a single molecules case.. please correct me! Anyway results for gas phase are close and this is not the issue now). Liquid phase: Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-27.3083 0.01 0.296591 -0.0389173 (kJ/mol) Coulomb (SR)6.00527 0.0074 0.122878 0.00576827 (kJ/mol) Coul. recip.5.59559 0.0032 0.0557413 0.00316957 (kJ/mol) Potential *34.6779 *0.0251.03468 -0.11177 (kJ/mol) Total Energy86.4044 0.0261.44353 -0.112587 (kJ/mol) *one single molecule in gas phase* Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-2.24473 0.073 1.292 0.342696 (kJ/mol) Coulomb (SR)11.5723 0.552.17577 -2.33224 (kJ/mol) Potential * 59.244 * 0.9410.97566.35631 (kJ/mol) Total Energy106.647 115.48286.78792 (kJ/mol) *equivalent number of molecules as in liquid* ( large box 20 nm) Statistics over 101 steps [ 0. through 2000. ps ], 4 data sets All statistics are over 11 points Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-2.16367 0.053 0.171542 0.374027 (kJ/mol) Coulomb (SR)11.2894 0.230.49105 -1.44437 (kJ/mol) Potential * 63.2369*1.12.472117.69756 (kJ/mol) Total Energy114.3371.12.655477.72258 (kJ/mol) Since pbc is set to NO molecules leave the box and I dont know if this all right. I hope the difference is acceptable...! For pbc = no there is no box. 0- I am going to do the same calculation but for some polymers solvated in the alkane. For binary system do I need to look at nonboded terms? and then run a simulation for a single polymer in vacuum? Can you please provide me with a recipe for Delta Hvap of the solute in a solvent? The method for calculating heat of vaporization is not dependent upon the contents of the system; it is a fundamental thermodynamic definition. Heat of vaporization is not something that can be calculated from a solute in a solvent. You can calculate DHvap for a particular system, but not some subset of that system. Thanks Justin. I am interested in the energy required to vaporize the solute in a particular solvent not the whole DHvap of the mixture. do you think this can be achieved by calculating nonbonded energies between solute and solvent? ( defining energy groups ..) 1- If I want to look at nonboded interactions only, do I have to add Coul. recip. to [ LJ (SR) + Coulomb (SR) ] ? The PME-related terms contain both solute-solvent, solvent-solvent, and potentially solute-solute terms (depending on the size and nature of the solute), so trying to interpret this term in some pairwise fashion is an exercise in futility. What you mean is when one uses PME interaction energies between components can not be decomposed? So the energy groups I defined to extract nonbonded energies are not giving correct values? Sofar I have been defining energy groups to calculate nonbonded terms between components _A-A A_B... I hope I have not been doing thing wrongly! Please help me out! Thanks, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http
[gmx-users] Heat of vap
Dear David, I followed your instructions and calculated Heat of vaporization of my alkane once with one molecule in gas phase (no cutoff) and once with equivalent number of molecules as in liquid phase as Justin suggested. Results are as follows: *one single molecule in gas phase* Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-2.24473 0.073 1.292 0.342696 (kJ/mol) Coulomb (SR)11.5723 0.552.17577 -2.33224 (kJ/mol) Potential * 59.244 * 0.9410.97566.35631 (kJ/mol) Total Energy106.647 115.48286.78792 (kJ/mol) *equivalent number of molecules as in liquid* ( large box 20 nm) Statistics over 101 steps [ 0. through 2000. ps ], 4 data sets All statistics are over 11 points Energy Average Err.Est. RMSD Tot-Drift --- LJ (SR)-2.16367 0.053 0.171542 0.374027 (kJ/mol) Coulomb (SR)11.2894 0.230.49105 -1.44437 (kJ/mol) Potential * 63.2369*1.12.472117.69756 (kJ/mol) Total Energy114.3371.12.655477.72258 (kJ/mol) Since pbc is set to NO molecules leave the box and I dont know if this all right. I hope the difference is acceptable...! 0- I am going to do the same calculation but for some polymers solvated in the alkane. For binary system do I need to look at nonboded terms? and then run a simulation for a single polymer in vacuum? Can you please provide me with a recipe for Delta Hvap of the solute in a solvent? 1- If I want to look at nonboded interactions only, do I have to add Coul. recip. to [ LJ (SR) + Coulomb (SR) ] ? coulombtype = PME vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0.9 rlist = 1.2 rcoulomb= 1.2 rvdw= 1.0 ff: OPLSAA Appreciate your help, Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Heat of vap
in your mail: On 30 March 2011 15:30, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: Dear all, I intend to obtain vaporization heat per volume for a /pure alkane system/. Here is the steps I am taking. Please correct me. 1- Obtain total energy of system (kinetic+potential) and divide by number of molecules to obtain energy per mol of molecules. g_energy -f *.edr -nmol XXX 2- Obtain total energy of a single molecule (use pbc). 3- Subtract step 2 from step 1. 4- Divide by simulation box volume. My questions is: in step 2 : what should be the box size? The same size as in 1 or it does not matter? (step 1 is done for the actual denstiy) More troubling, how does one define the energy of a molecule? If you use any sort of long-range algorithms (especially PME, but also dispersion correction), you can't simply decompose the system like this. Thanks Justin and David. I have been trying to find the article in which this has been presented. If you have time Please see page 5937, right column, equation 11. I think I made a mistake and I dont have to include kinetic energy, Only nonboded energies!? http://pubs.acs.org/doi/pdfplus/10.1021/jp0707539 In the derivation of recent Gromos96 parameter sets, the heat of vaporization is quite simple: DHvap = Ugas - Uliq + RT 1- So Uliq is the total energy or only potential (no kinetic) 2- How can I compute Ugas? I have liquid now... Thank you, Regards, -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Heat of vap
On 31 March 2011 12:58, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: in your mail: On 30 March 2011 15:30, Justin A. Lemkul jalem...@vt.edu mailto: jalem...@vt.edu wrote: Elisabeth wrote: Dear all, I intend to obtain vaporization heat per volume for a /pure alkane system/. Here is the steps I am taking. Please correct me. 1- Obtain total energy of system (kinetic+potential) and divide by number of molecules to obtain energy per mol of molecules. g_energy -f *.edr -nmol XXX 2- Obtain total energy of a single molecule (use pbc). 3- Subtract step 2 from step 1. 4- Divide by simulation box volume. My questions is: in step 2 : what should be the box size? The same size as in 1 or it does not matter? (step 1 is done for the actual denstiy) More troubling, how does one define the energy of a molecule? If you use any sort of long-range algorithms (especially PME, but also dispersion correction), you can't simply decompose the system like this. Thanks Justin and David. I have been trying to find the article in which this has been presented. If you have time Please see page 5937, right column, equation 11. I think I made a mistake and I dont have to include kinetic energy, Only nonboded energies!? http://pubs.acs.org/doi/pdfplus/10.1021/jp0707539 What is cohesive energy and how does it relate to the quantity you're trying to calculate? It is delta Hvap/volume. It is directly related to Hvap. What is happening is that they are calculating nonbonded energy of some chains, divide by number of chains and substract from nonbonded energy of a single chain in vacuum. These are the steps I wrote in my first post but I think I should not have included kinetic and should just look at LJ-SR and Coulomb-SR. I am using PME..If I remember correctly LR is included in Coulomb-SR and can not get decomposed? But I dont think this doesnt matter since if I am to take nonbonded energies this should not hurt,,, Please comment ... In the derivation of recent Gromos96 parameter sets, the heat of vaporization is quite simple: DHvap = Ugas - Uliq + RT 1- So Uliq is the total energy or only potential (no kinetic) Potential. 2- How can I compute Ugas? I have liquid now... Run a simulation in the gas phase. Sorry, but how can I do this? :( I have box of molecules with density of actual liquid..How can I shift to gas phase ..I mean how many molecules I need to keep in the box.. Many thanks... -Justin Thank you, Regards, -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu http://vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] index file
Dear gromacs users, I am planning to produce C-C rdf for my system and I have difficulty making the index file with a group of all C in the system, [Carbon] . .If there are thousands of atoms in the system do I have to go through the top file and note down C atom numbers and use those numbers to make a group in the index file? What is the easiest way of doing this using make_ndx features.. Thank you for your hints.. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Thermal expansion coeff.
Dear all, I asked this question a few days ago but did not receive a solid answer. Has anyone tried to produce the properties below? Can you comment on accuracy of these properties using g_energy? Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp Isothermal compressibility: Vol, Temp Adiabatic bulk modulus: Vol, Temp Thanks, Kind regards, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Heat of vap
Dear all, I intend to obtain vaporization heat per volume for a *pure alkane system*. Here is the steps I am taking. Please correct me. 1- Obtain total energy of system (kinetic+potential) and divide by number of molecules to obtain energy per mol of molecules. g_energy -f *.edr -nmol XXX 2- Obtain total energy of a single molecule (use pbc). 3- Subtract step 2 from step 1. 4- Divide by simulation box volume. My questions is: in step 2 : what should be the box size? The same size as in 1 or it does not matter? (step 1 is done for the actual denstiy) Thanking you all, Regards, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] heat capacity in 4.5.3 and 4.0.7
Thanks David for replying me. I appreciate your attention to comment on the following short inquiries. 1- g_energy says: you will need to give the number of constraints per molecule –nconstr and for water this is 3. Can you elaborate on how to set this option. When constraints= all-bonds how to set -nconstr for a simlpe alkane for instance? It refers to number of atoms of each molecule on which constraints are applied? 2- Aside from heat capacities, does g_energy report the following properties correctly? Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp Isothermal compressibility: Vol, Temp Adiabatic bulk modulus: Vol, Temp 3- Since volume appears for calculation of the above properties , does this imply ONLY NPT results in these properties? (because in NPT g_energy reports volume). 4- For thermal expansion: I read one of your earlier messages saying that enthalpy = Etot(total energy) +pv, [where Etot is potential+kinetic (in g_energy list)] and for pv , take ref_p and not measure p (that g_energy reports). If the system is equilibrated to ref_p, then why not taking pv from g_energy? you mean product of pv thatg_energy gives is not correct? Many thanks for you time! Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] heat capacity in 4.5.3 and 4.0.7
Dear all, I am trying to calculate Cv of a pure alkane 125 molecules with actual density, and here is what I have done: 1- Cv: Run NVT g_energyselect T and Etot=total energy As I issue g_energy *.edr Heat Capacity Cv: 12.4787 J/mol K (factor = 0.000369503) And when I choose –nmol XXX Heat Capacity Cv: 13.4007 J/mol K (factor = 0.0462137) - Do I have to use –nmol? What would be the unit with –nmol ? and as I issue g_energy using 4.5.3 Temperature dependent fluctuation properties at T = 299.987. #constr/mol = 0 Heat capacity at constant volume Cv:341.734 J/mol K and with -nconstr # of bonds Temperature dependent fluctuation properties at T = 299.987. #constr/mol = 19 Heat capacity at constant volume Cv:262.746 J/mol K - From this I am guessing what 4.0.7 reports is not correct!? but still what I am getting from 4.5.3 is far from actual reported value (around 170). I guess I am missing something..! or Is the difference acceptable? g_energy says: you will need to give the number of constraints per molecule –nconstr - so I need to count how many bonds are fixed in the molecule as I am setting constraints= all-bonds? The way to this is counting no. of bonds? What about macromolecules? Is there an easier way of finding number of constraints? 2- Cp: I did NPT, g_energy *.edr -nmol XXX -nconstr (4.5.3) I want to select T and enthalpy but since enthalpy is not on the list I select total-energy (which is Etot I think) and pv terms but g_energy is still giving Cv! PropertyEnergy terms needed --- Heat capacity Cp (NPT sims):Enthalpy, Temp Energy Average Err.Est. RMSD Tot-Drift --- Total Energy86.1571 0.0851.55565 -0.0944531 (kJ/mol) Temperature 300.074 0.075.77219 0.318146 (kJ/mol) pV 15.8128 0.0611080.82 -0.249203 (kJ/mol) Temperature dependent fluctuation properties at T = 300.074. #constr/mol = 19 Heat capacity at constant volume Cv:325.067 J/mol K Can you please guide me where the problem lies? --- again Cv : to compare with constraints=none but still -nconstr should be on otherwise I values are more off. Total Energy132.418 0.0511.77314 0.144794 (kJ/mol) Temperature 300.008 0.0524.502090.15866 (kJ/mol) Temperature dependent fluctuation properties at T = 300.008. #constr/mol = 20 Heat capacity at constant volume Cv:442.016 J/mol K 3- Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp Isothermal compressibility: Vol, Temp Adiabatic bulk modulus: Vol, Temp Since volume appears for calculation of the above properties , does this imply only NPT results in these properties? and wanted to make sure what 4.5.3 reports is reliable.. Many thanks for you time! Best, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] surface tension vs. system size
Dear all, I performed surface tension calculations vs. different system size (NVT, mdp file is included at the end). The reported surface tension for hydrocarbon I am studying is 18 mN/m at 20C. I am getting ~ 175 bar nm from g_energy which means surface tension of around 9 mN/m. For box size 3 3 6 nm I get 202 bar nm which is 10 mN/m and is closer to the reported one. 0-Are my results reasonable? Is it possible to get closer to the actual surf. ten? or the current difference can not be improved by molecular dynamics? 1- I am just wondering how I can decide on the system size given runs performed (shown below). Which system I can stick to? 2- Can I conclude that surface tension I am getting is equilibrated one? RMSD is much larger than average!! There are 125 molecules in a 3 3 3 nm box initially. *F*or 6 6 ? runs: 3 3 3 was replicated in X and Y. (that is 4*125 molecules) *box size 6 6 18 nm (molecules fill up volume of 6 6 3 and Z direction is increased to 18 so total size is 6 6 18) * Statistics over 898601 steps [ *0. through 1797.2001 ps ]*, 1 data sets All statistics are over 898601 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen175.322.91813.06 -17.41 (kJ/mol) Statistics over 818001 steps [ *0. through 1636.0001 ps* ], 1 data sets All statistics are over 818001 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 175.4824.91814.64 -19.6398 (kJ/mol) Statistics over 569501 steps [ *500. through 1639.0001 ps* ], 1 data sets All statistics are over 569501 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 175.459 51815.85 -31.5002 (kJ/mol) Statistics over 321001 steps [* 1000.0001 through 1642.0001 ps ]*, 1 data sets All statistics are over 321001 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 168.1283.91812.23 -12.0096 (kJ/mol) *box size 6 6 8 nm. (molecules fill up volume of 6 6 3 and Z direction is increased to 8 so total size is 6 6 8) *Statistics over 101 steps [* 0. through 2000.0001 ps *], 1 data sets All statistics are over 101 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 177.2743.91819.27 -0.465725 (kJ/mol)* * Statistics over 750001 steps *[ 500. through 2000.0001 ps ],* 1 data sets All statistics are over 750001 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 176.2473.71817.8812.2119 (kJ/mol) * * - *box size 6 6 6 nm. (molecules fill up volume of 6 6 3 and Z direction is increased to 6 so total size is 6 6 6) *Statistics over 532601 steps [ 0. through 1065.2001 ps ], 1 data sets All statistics are over 532601 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 175.1372.21816.63 -2.51592 (kJ/mol) Statistics over 282601 steps [ *500. through 1065.2001 ps* ], 1 data sets All statistics are over 282601 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 178.4137.61815.25 -21.8612 (kJ/mol) - *box size 3 3 9 nm. (molecules fill up volume of 3 3 3 and Z direction is increased to 9 so total size is 3 3 9)* Statistics over 101 steps [ 0. through 2000.0001 ps ], 1 data sets All statistics are over 101 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 180.9858.23621.5518.4731 (kJ/mol) Statistics over 750001 steps [ 500. through 2000.0001 ps ], 1 data sets All statistics are over 750001
Re: [gmx-users] unit of surface tension, KJ/mol?
Thanks for your help. :) On 23 March 2011 07:47, Justin A. Lemkul jalem...@vt.edu wrote: Mark Abraham wrote: On 23/03/2011 3:39 PM, Elisabeth wrote: Thanks Mark for the hints. I am using a pre equilibrated box for both runs. Can you again explain how Can I interpret this from fluctuations from two runs? Below is the resutls from two versions: Thanks! IIRC the contents of the .edr file didn't change across these versions, but the reporting of computed quantities did, so I'd make my life simpler and use only the 4.5.3 g_energy. Actually, 4.5.3 may have the bug: http://redmine.gromacs.org/issues/696 The quantities in the .edr file are correct, but the output from g_energy is nonsense. It was fixed for 4.5.4. -Justin Both simulations have a much higher value for the root-mean-squared deviation from the average (i.e. standard deviation) than for the average. So that means your data set looks something like 1 20 6543 200 23444 23434 15 500 3444. Look at the .xvg file that is produced. You need a heck of a lot of such data to have confidence that your sample average reflects the true average. If you haven't got that pile of data, then you haven't observed the average with any confidence. This is normal for quantities computed from fluctuations in atomic positions. They're macroscopic quantities, and have to be observed over a lot of microscopic configurations. Be sure to check visually that your system is doing what you hope it is doing. Mark 4.0.7 Statistics over 733801 steps [ 0. thru 1467.6001 ps ], 1 data sets All averages are exact over 733801 steps Energy Average RMSD Fluct. Drift Tot-Drift --- #Surf*SurfTen 203.7783623.993623.98 -0.00203751 -2.99026 Statistics over 483401 steps [ 500. thru 1466.8000 ps ], 1 data sets All averages are exact over 483401 steps Energy Average RMSD Fluct. Drift Tot-Drift --- #Surf*SurfTen 203.1293612.583612.28 0.168686 163.086 - 4.5.3 Statistics over 732501 steps [ 1. through 1466. ps ], 1 data sets All statistics are over 73251 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 2374.4957021081.3 -3251.38 (bar nm) Statistics over 384501 steps [ 1000. through 1769. ps ], 1 data sets All statistics are over 38451 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 2610.25 100021109.94638.38 (bar nm) Statistics over 483001 steps [ 500. through 1466. ps ], 1 data sets All statistics are over 48301 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 1645.0739020288.3 -174.289 (bar nm) On 22 March 2011 23:56, Mark Abraham mark.abra...@anu.edu.au mailto: mark.abra...@anu.edu.au wrote: On 23/03/11, *Elisabeth * katesed...@gmail.com mailto:katesed...@gmail.com wrote: On 22 March 2011 22:46, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Elisabeth wrote: On 22 March 2011 22:31, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu mailto:jalem...@vt.edu mailto:jalem...@vt.edu wrote: Elisabeth wrote: Hello, I did two simulations on the same system using versions 4.0.7 and 4.5.3. It seems like the unit of surface tension is not the same in these versions because I am getting ~250 KJ/mol an in 4.0.7 and ~ 5000bar nm in 4.5.3?! How KJ/mol can be converted into bar nm? Can anyone help please. Is this from g_energy output? In past versions, everything was printed as kJ/mol, even quantities that obviously weren't, like temperature, pressure, etc. Yes, so why are the results so different. I am using exactly the same mdp file.! Any pressure-related quantity is going to be subject to enormous fluctuations. This has been discussed within the last few days. Without seeing the .mdp file and a description of the system, it's hard for anyone
Re: [gmx-users] unit of surface tension, KJ/mol?
On 22 March 2011 22:31, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: Hello, I did two simulations on the same system using versions 4.0.7 and 4.5.3. It seems like the unit of surface tension is not the same in these versions because I am getting ~250 KJ/mol an in 4.0.7 and ~ 5000bar nm in 4.5.3?! How KJ/mol can be converted into bar nm? Can anyone help please. Is this from g_energy output? In past versions, everything was printed as kJ/mol, even quantities that obviously weren't, like temperature, pressure, etc. Yes, so why are the results so different. I am using exactly the same mdp file.! -Justin Thanks -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] unit of surface tension, KJ/mol?
On 22 March 2011 22:46, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: On 22 March 2011 22:31, Justin A. Lemkul jalem...@vt.edu mailto: jalem...@vt.edu wrote: Elisabeth wrote: Hello, I did two simulations on the same system using versions 4.0.7 and 4.5.3. It seems like the unit of surface tension is not the same in these versions because I am getting ~250 KJ/mol an in 4.0.7 and ~ 5000bar nm in 4.5.3?! How KJ/mol can be converted into bar nm? Can anyone help please. Is this from g_energy output? In past versions, everything was printed as kJ/mol, even quantities that obviously weren't, like temperature, pressure, etc. Yes, so why are the results so different. I am using exactly the same mdp file.! Any pressure-related quantity is going to be subject to enormous fluctuations. This has been discussed within the last few days. Without seeing the .mdp file and a description of the system, it's hard for anyone to comment on what the results might mean. Thanks. I am working on a pure alkane system, in a box of 3X3X3 which is extended in Z to create liq/air interface. That is 3X3X6 . What could be reason for such a big difference in results from two version? Thanks alot! integrator = md dt = 0.002 nsteps = 100 nstenergy = 100 nstxout = 100 nstlist = 10 ns_type = grid coulombtype = PME vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0 rlist = 1.1 rcoulomb= 1.1 rvdw= 1.0 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 ;optimize_fft = yes Tcoupl = v-rescale tc-grps = System tau_t = 0.1 ref_t = 300 pbc = xyz gen_vel = yes gen_temp= 300.0 gen_seed= 173529 constraints = all-bonds constraint-algorithm = lincs -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] unit of surface tension, KJ/mol?
Thanks Mark for the hints. I am using a pre equilibrated box for both runs. Can you again explain how Can I interpret this from fluctuations from two runs? Below is the resutls from two versions: Thanks! 4.0.7 Statistics over 733801 steps [ 0. thru 1467.6001 ps ], 1 data sets All averages are exact over 733801 steps Energy Average RMSD Fluct. Drift Tot-Drift --- #Surf*SurfTen 203.7783623.993623.98 -0.00203751 -2.99026 Statistics over 483401 steps [ 500. thru 1466.8000 ps ], 1 data sets All averages are exact over 483401 steps Energy Average RMSD Fluct. Drift Tot-Drift --- #Surf*SurfTen 203.1293612.583612.28 0.168686 163.086 - 4.5.3 Statistics over 732501 steps [ 1. through 1466. ps ], 1 data sets All statistics are over 73251 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 2374.4957021081.3 -3251.38 (bar nm) Statistics over 384501 steps [ 1000. through 1769. ps ], 1 data sets All statistics are over 38451 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 2610.25 100021109.94638.38 (bar nm) Statistics over 483001 steps [ 500. through 1466. ps ], 1 data sets All statistics are over 48301 points Energy Average Err.Est. RMSD Tot-Drift --- #Surf*SurfTen 1645.0739020288.3 -174.289 (bar nm) On 22 March 2011 23:56, Mark Abraham mark.abra...@anu.edu.au wrote: On 23/03/11, *Elisabeth * katesed...@gmail.com wrote: On 22 March 2011 22:46, Justin A. Lemkul jalem...@vt.edu wrote: Elisabeth wrote: On 22 March 2011 22:31, Justin A. Lemkul jalem...@vt.edu mailto: jalem...@vt.edu wrote: Elisabeth wrote: Hello, I did two simulations on the same system using versions 4.0.7 and 4.5.3. It seems like the unit of surface tension is not the same in these versions because I am getting ~250 KJ/mol an in 4.0.7 and ~ 5000bar nm in 4.5.3?! How KJ/mol can be converted into bar nm? Can anyone help please. Is this from g_energy output? In past versions, everything was printed as kJ/mol, even quantities that obviously weren't, like temperature, pressure, etc. Yes, so why are the results so different. I am using exactly the same mdp file.! Any pressure-related quantity is going to be subject to enormous fluctuations. This has been discussed within the last few days. Without seeing the .mdp file and a description of the system, it's hard for anyone to comment on what the results might mean. Thanks. I am working on a pure alkane system, in a box of 3X3X3 which is extended in Z to create liq/air interface. That is 3X3X6 . What could be reason for such a big difference in results from two version? Thanks alot! Look at the size of the fluctuations printed by g_energy. If they're comparable or larger than your differences, then you are not observing a statistically significant difference. Over just 2ns of a small system, that's likely to be the case. Damping such fluctuations requires large simulations or long times or both. You are also generating velocities at the start of this run, so probably you are including spurious measurements on non-equilibrated configurations in your statistics. If you'd thought to tell us your command lines, I wouldn't be guessing :-) You can fix this with g_energy -b, but really you should separate your equilibration from your production run, even if you don't really need to, so that all your workflows have similar properties. Then when you come back to repeat some analysis in a month's time, you don't have to remember to ignore the first 1ns of this simulation... Mark integrator = md dt = 0.002 nsteps = 100 nstenergy = 100 nstxout = 100 nstlist = 10 ns_type = grid coulombtype = PME vdw-type= Shift rcoulomb-switch = 0 rvdw-switch = 0 rlist = 1.1 rcoulomb= 1.1 rvdw= 1.0 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 ;optimize_fft
[gmx-users] surface tension
Hello, Thank you for your answer. 1- If I am right I have to increase the length in two directions rather than one, to create a plane parallel to XY for example? 2- Can you please give me an idea on how many molecules I need to have in the box and also what should be the thickness of layer? I have now 3nm X 9 X 9 dimensions. That is thickness of 3nm. What I did was replicating a 3nm box using genconf -nbox 3 3 1. I dont know what is the correct way of creating a layer for surface tension calculation. I appreciate any comments about number of molecules, box dimensions for such a study. 3- my last question is how can I make sure surface tension reported by g_energy is the equilibrated one. RMSD is very big compared to surf. ten. ! Thanks for your time. Elisabeth ** if you are interested in the surface tension of a pure liquid, which I assume is true from your message, then you need to create at least one surface, since periodic boundary conditions make the model system infinite, i.e., without a surface whatsoever. the easiest way to make that happen is to increase the length of the box in one direction, say the z direction. that way you will end up with a system that resemble a (thin) liquid film with vacuum below and above, meaning that you now have two surfaces. run a regular simulation (NVT) e use g_energy to get the surface tension. btw: as any other pressure related property, fluctuations are huge. best Andre On Wed, Mar 9, 2011 at 12:25 PM, Elisabeth katesed...@gmail.com wrote: Dear gmx users, Since I am new to surface tension topic I need to ask very trivial questions. Please help me out with these simple questions. As a starting point I am going to calculate surface tension of a pure alkane in a cubic box and compare with experimental values. 1- g_energy is giving #Surf*SurfTen by default. On the other hand surface tension can be obtained by gamma = (Pzz - (Pxx+Pyy)/2) / Lz. i.e Pres-XX-(bar), Pres-YY(bar), Pres--(bar) Can anyone tell me what the difference between these two is? 2- In pressure coupling settings there is surface_tension option which I guess is applicable where surface tension needs to be kept fixed. If one want to calculate surface tension I dont think this option make sense. Am I right? 3- I am using the following setting: I calculate the average for a 2ns run and different start times as shown below. Although T, P and other quantities are equilibrated after 200ps, surface tension is not giving a constant value. Is that because I am not using berenden P coupling? (As mentioned in the manual surface tension works with berendsen) Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 1 compressibility = 4.5e-5 ref_p = 40 time period for which average is calculated Average RMSD Fluct. Drift Tot-Drift --- 1-2000 ps run: #Surf*SurfTen6.438443588.74 3588.35 0.091406182.721 500-2000 ps#Surf*SurfTen 12.85183605.72 3605.26 0.132126198.189 1000-2000ps #Surf*SurfTen18.88213610.97 3610.80.11819118.191 1500-2000ps #Surf*SurfTen 23.00723585.51 3584.93 -0.444037 -222.019 4- Assuming I am getting surface tension for a cubic box, to compare this with reported values in literature I need to divide by 6 (no. of surfaces)? 5- Does box six affect the results? (mine is 3.3 nm ). Thank you, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: surface tension
Dear Andre, Thanks for the helpful information. I need to do some text reading to understand the periodic BC effect you are talking about. I dont see why increasing length in z direction does not lead to periodic BC in z and only for x, y ? does that mean the thickness of layer would be the Z dimension? then how much increase in one direction is reasonable? (If I have a 2 nm box). Also Can you please introduce some text book? Thank you, Best regards, 2011/3/15 André Farias de Moura mo...@ufscar.br Dear Elisabeth, actually, it is the other way around, you need increase the box length in one direction, thus keeping periodic boundary conditions in the other two directions while a (infinitely periodic) surface is created. and notice that using genconf with -nbox 3 3 1 will increase your system but will not make it a surface system, unless you increase the box length in one direction. as regards the size, the larger the model system, the smaller should be the fluctuations. but mind that you should increase the size by a few order of magnitude for any noticeable decrease on the (huge) RMSD values you are getting. if you want to check the results for convergence maybe you could try either a block averaging or a running average (grace can do that for you). best regards Andre On Mon, Mar 14, 2011 at 7:57 PM, Elisabeth katesed...@gmail.com wrote: Hello, Thank you for your answer. 1- If I am right I have to increase the length in two directions rather than one, to create a plane parallel to XY for example? 2- Can you please give me an idea on how many molecules I need to have in the box and also what should be the thickness of layer? I have now 3nm X 9 X 9 dimensions. That is thickness of 3nm. What I did was replicating a 3nm box using genconf -nbox 3 3 1. I dont know what is the correct way of creating a layer for surface tension calculation. I appreciate any comments about number of molecules, box dimensions for such a study. 3- my last question is how can I make sure surface tension reported by g_energy is the equilibrated one. RMSD is very big compared to surf. ten. ! Thanks for your time. Elisabeth ** if you are interested in the surface tension of a pure liquid, which I assume is true from your message, then you need to create at least one surface, since periodic boundary conditions make the model system infinite, i.e., without a surface whatsoever. the easiest way to make that happen is to increase the length of the box in one direction, say the z direction. that way you will end up with a system that resemble a (thin) liquid film with vacuum below and above, meaning that you now have two surfaces. run a regular simulation (NVT) e use g_energy to get the surface tension. btw: as any other pressure related property, fluctuations are huge. best Andre On Wed, Mar 9, 2011 at 12:25 PM, Elisabeth katesed...@gmail.com wrote: Dear gmx users, Since I am new to surface tension topic I need to ask very trivial questions. Please help me out with these simple questions. As a starting point I am going to calculate surface tension of a pure alkane in a cubic box and compare with experimental values. 1- g_energy is giving #Surf*SurfTen by default. On the other hand surface tension can be obtained by gamma = (Pzz - (Pxx+Pyy)/2) / Lz. i.e Pres-XX-(bar), Pres-YY(bar), Pres--(bar) Can anyone tell me what the difference between these two is? 2- In pressure coupling settings there is surface_tension option which I guess is applicable where surface tension needs to be kept fixed. If one want to calculate surface tension I dont think this option make sense. Am I right? 3- I am using the following setting: I calculate the average for a 2ns run and different start times as shown below. Although T, P and other quantities are equilibrated after 200ps, surface tension is not giving a constant value. Is that because I am not using berenden P coupling? (As mentioned in the manual surface tension works with berendsen) Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 1 compressibility = 4.5e-5 ref_p = 40 time period for which average is calculated Average RMSD Fluct. Drift Tot-Drift --- 1-2000 ps run: #Surf*SurfTen6.438443588.74 3588.35 0.091406182.721 500-2000 ps#Surf*SurfTen 12.85183605.72 3605.26 0.132126198.189 1000-2000ps #Surf*SurfTen18.88213610.97 3610.80.11819118.191 1500-2000ps #Surf*SurfTen 23.00723585.51 3584.93 -0.444037 -222.019 4- Assuming I am getting surface tension
[gmx-users] surface tension in gmx
Dear gmx users, Since I am new to surface tension topic I need to ask very trivial questions. Please help me out with these simple questions. As a starting point I am going to calculate surface tension of a pure alkane in a cubic box and compare with experimental values. 1- g_energy is giving #Surf*SurfTen by default. On the other hand surface tension can be obtained by gamma = (Pzz - (Pxx+Pyy)/2) / Lz. i.e Pres-XX-(bar), Pres-YY(bar), Pres--(bar) Can anyone tell me what the difference between these two is? 2- In pressure coupling settings there is surface_tension option which I guess is applicable where surface tension needs to be kept fixed. If one want to calculate surface tension I dont think this option make sense. Am I right? 3- I am using the following setting: I calculate the average for a 2ns run and different start times as shown below. Although T, P and other quantities are equilibrated after 200ps, surface tension is not giving a constant value. Is that because I am not using berenden P coupling? (As mentioned in the manual surface tension works with berendsen) Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p = 1 1 compressibility = 4.5e-5 ref_p = 40 time period for which average is calculated Average RMSD Fluct. Drift Tot-Drift --- 1-2000 ps run: #Surf*SurfTen6.438443588.74 3588.35 0.091406182.721 500-2000 ps#Surf*SurfTen 12.85183605.72 3605.26 0.132126198.189 1000-2000ps #Surf*SurfTen18.88213610.97 3610.80.11819118.191 1500-2000ps #Surf*SurfTen 23.00723585.51 3584.93 -0.444037 -222.019 4- Assuming I am getting surface tension for a cubic box, to compare this with reported values in literature I need to divide by 6 (no. of surfaces)? 5- Does box six affect the results? (mine is 3.3 nm ). Thank you, -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
