[gmx-users] Fwd: PE compressing - for reference only
Vitaly Chaban
vvchaban at gmail.com
Sat Jul 24 15:37:31 CEST 2010
Dear Moeed:
Please find attached the simulation results as well as the input files
for your PE system. I started the simulation from your 8PE
configuration and compressed it to the density of ~890 kg/m3. It is
rather good since ranges between the experimental values -
http://en.wikipedia.org/wiki/Polyethylene .The system is also in
equilibrium (please see the attached plot). I think this is what you
need to start productive runs.
You might also want to adjust charges but this is up to you. I suppose
the average density can become slightly larger if you increase the
system size because now the system is very modest as compared to the
size of one PE particle.
> 1- In general in most of the NPT runs I am not reaching the reference
> pressure ( 50 ,30 ..bar). I get very low values around zero or 1.
I've just performed the simulation at 30 bar during 2 000 ps and got
Energy Average RMSD Fluct. Drift Tot-Drift
-------------------------------------------------------------------------------
Pressure (bar) 30.1703 1373.51 1368.76 -0.790566
-395.285
Thus, no problems with NPT... I also think you do not NVT at all in this study.
> 2-About bond constraints: Actually, in literature I see it is common to fix
> the bond distance. Also in one of the GROMACS tutorials I have seen the
> reasoning that" since bond vibrations are of quantum nature they can not be
> treated classically". The reason you mentioned that I am noy using a good
> force filed is because of using bond constraints or ...? Do I need to remove
> constraints?
The fixation of the interatomic distances is not so common. The
general rule is not to use constraints&restraints if you can avoid
their usage. In my opinion, your system need not constaints. I did not
just used them on your system (see MDP attached).
The force field looked suspicious since you got LINCS warning. With a
correct setup and interaction parameters, no warnings should be
present.
> 3- I have another inquiry about treating electrostatics contribution in my
> system. My systems as you have noticed involves only C and H. Later I will
> have to solvate this PE chain in HExane which is again apolar (only partial
> charges). It is not still clear to me whether I need to include
> electrostatic energy terms or I need to consider only vdw for calculating
> interaction energies
>From the scientific point of view, I'd suggest to run both cases and
then conclude about the role of electrostatics in your particular
polymer system. Generally, you need electrostatics in order to
represent all the non-bonded interactions correctly, even if the
partial charges are negligible. I do not like SHIFT for
electrostatics, please use PME as I recently did with your system
(attached MDP).
> 4-the value of compressibility I am using in mdp file is 4.5e-5 which is for
> water. I have not found this valued for my polymer yet.
This value is not critical. Usually, the guys use the same value
(almost) for all systems. Read about Parrinello-Rahman pressure
coupling for the more comprehensive insight on how to get
compressibility.
Also, please pay attention to the modified MDP file. These parameters
look better for your particular system.
For the productive runs, I'd suggest you to enlarge the simulation
box, so it is significantly longer or at least not shorter than the PE
chain. I also suppose it will give some increase of density since the
voids due to polymer complicated structure will become smaller as
compared to the overall box volume.
Good luck!
Dr. Vitaly Chaban
On Sat, Jul 24, 2010 at 5:42 AM, Moeed <lecielll at googlemail.com> wrote:
> Dear Dr. Chaban,
>
> First of all, I thank you for you kind attention. I am a PhD student at
> McGill and my topic revolves around computational thermodynamics of polymer
> solutions. I appreciate any comment or suggestion...
>
> ( I am focusing on calculating interaction energies between PE HExane and
> using these values in thermodynamic model. At the moment I have only PE and
> would like to study a pure polymer system. I would like also to look into
> free energy of solvation, rms, Rg, .. and very importantly surface tension
> calculations)
>
> I noticed that in either NPT or NVT trials for a long chain (PE60.gro)
> equilibrium state corresponds to more extended conformations. I noticed the
> only way to approach the volume I am after is forget the equilibrium state
> (as if simulation is long chain tends to extend) and I looked for the
> globule-like structure so the volume it occupies is minimal. For instance I
> did a NPT for 100ps and viewing the trajectory I noticed at 30 ps molecule
> has the convoluted structure so in the next trial I picked simul. time of 30
> ps to capture the structure I want (PE60-P150-30ps, pressure 150) , then
> used editconf to reduce the size (editconf888 8*8*8nm) energy minimized the
> structure and performed the next NPT with this initial structure. Doing so,
> I could approach the box size of 2.26 nm (desired size is 1.8nm).
>
> File names I have picked make it easy to follow the steps (especially if you
> view the files using "List option" or "details option" in your directory) I
> have taken from the initial PE chain (PE60.gro)
>
> PE60-P150-30ps-editconf888-P50-100ps-NVT140ps-editconf444center_bfpr.gro
> this is still a single chain in box of size 4 4 4 and _b4pr indicates the
> last step has been EM.
>
> after this I did two more NPT runs:
>
> PE60-P150-30ps-editconf888-P50-100ps-NVT140ps-editconf444center_bfpr-NPT120-100ps-NPT50-100p.gro
>> this is renamed to PE60-box2.26
>
> and got box size of 2.26nm
>
> Since it was not possible to compress the system further I replicated this
> molecule and energy minimized the replicated system to 8 chains
> (PE60-box2.26-nbox222). Energy minimization gives reasonable force but
> potential is about 3000. I dont know if this is reasonable for such a packed
> system (PE60-box2.26-nbox222_bfpr)
>
> Then employed NPT to compress again to a volume a bit smaller (3.34nm) than
> what I need and used editconf to attain the exact density I am after( which
> is size of 3.6nm after replicating 8 chains). Below is the output file of
> EM of the replicated and compressed system:
>
> PE60-box2.26-nbox222_bfpr-P50-180ps-editconf3.6c_b4pr
>
> I am getting some LINC warnings but at the end force and potential energies
> seem to be OK. PLease let me know if anything is going wrong.
> :
> Steepest Descents:
> Tolerance (Fmax) = 1.00000e+03
> Number of steps = 400
> Step= 0, Dmax= 1.0e-02 nm, Epot= 1.32635e+08 Fmax= 2.29301e+10, atom=
> 1004
> Step= 1, Dmax= 1.0e-02 nm, Epot= 2.22267e+07 Fmax= 1.39980e+09, atom=
> 1766
> Step= 2, Dmax= 1.2e-02 nm, Epot= 8.21648e+06 Fmax= 2.00629e+08, atom=
> 968
> Step= 3, Dmax= 1.4e-02 nm, Epot= 1.67342e+06 Fmax= 2.38551e+07, atom=
> 968
> Step= 4, Dmax= 1.7e-02 nm, Epot= 4.59581e+05 Fmax= 3.28726e+06, atom=
> 121
> Step= 5, Dmax= 2.1e-02 nm, Epot= 1.51948e+05 Fmax= 4.17261e+05, atom=
> 1766
> Step= 6, Dmax= 2.5e-02 nm, Epot= 4.95460e+04 Fmax= 1.64740e+05, atom=
> 121
> Step= 7, Dmax= 3.0e-02 nm, Epot= 3.20548e+04 Fmax= 8.22273e+04, atom=
> 968
> Step= 8, Dmax= 3.6e-02 nm, Epot= 2.16962e+04 Fmax= 1.86463e+04, atom=
> 961
> Step= 9, Dmax= 4.3e-02 nm, Epot= 1.29190e+04 Fmax= 3.92749e+03, atom=
> 1766
>
> Step 10, time 0.02 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.001238, max 0.029289 (between atoms 113 and 115)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> 113 114 33.6 0.1090 0.1103 0.1090
> 1765 1766 36.2 0.1090 0.1075 0.1090
> Step= 10, Dmax= 5.2e-02 nm, Epot= 8.04986e+03 Fmax= 3.49262e+04, atom=
> 826
>
> .
> .
> .Step 15, time 0.03 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.042648, max 1.350503 (between atoms 825 and 826)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> 122 124 38.1 0.1093 0.1174 0.1090
> 122 123 34.0 0.1092 0.1165 0.1090
> 119 121 41.9 0.1092 0.1253 0.1090
> 119 120 44.7 0.1091 0.1219 0.1090
> 116 119 33.7 0.1541 0.1512 0.1529
> 113 116 40.5 0.1540 0.1564 0.1529
> 113 115 55.8 0.1095 0.1294 0.1090
> 113 114 49.4 0.1093 0.1263 0.1090
> 110 112 50.6 0.1094 0.1209 0.1090
> 110 111 53.4 0.1094 0.1218 0.1090
> 107 110 39.4 0.1534 0.1666 0.1529
> 107 108 35.4 0.1090 0.1178 0.1090
> 966 969 33.5 0.1534 0.1669 0.1529
> 966 968 44.8 0.1090 0.1172 0.1090
> 966 967 66.8 0.1091 0.1221 0.1090
> 963 965 38.5 0.1091 0.1181 0.1090
> 963 964 37.1 0.1091 0.1179 0.1090
> 831 834 31.6 0.1538 0.1448 0.1529
> 831 833 37.3 0.1099 0.0977 0.1090
> 828 831 55.2 0.1553 0.1785 0.1529
> 828 829 84.4 0.1113 0.1348 0.1090
> 825 828 56.5 0.1555 0.2214 0.1529
> 825 827 89.3 0.1102 0.1823 0.1090
> 825 826 89.2 0.1085 0.2562 0.1090
> 822 824 65.3 0.1092 0.1602 0.1090
> 822 823 51.3 0.1092 0.1482 0.1090
> 819 822 44.6 0.1534 0.1902 0.1529
> 1768 1770 45.1 0.1107 0.1528 0.1090
> 1768 1769 38.5 0.1101 0.1401 0.1090
> 1765 1768 55.7 0.1510 0.1687 0.1529
> 1765 1767 88.7 0.1145 0.1453 0.1090
> 1765 1766 90.4 0.1128 0.1744 0.1090
> 1762 1765 36.3 0.1601 0.1635 0.1529
> 1762 1764 42.7 0.1120 0.1163 0.1090
> 1762 1763 48.8 0.1101 0.1144 0.1090
> 1759 1760 38.9 0.1102 0.1097 0.1090
> 1591 1594 33.0 0.1533 0.1547 0.1529
> 1591 1593 36.9 0.1092 0.1158 0.1090
> 1591 1592 36.1 0.1091 0.1153 0.1090
> 1588 1590 38.0 0.1092 0.1140 0.1090
> 1588 1589 39.0 0.1092 0.1141 0.1090
> 1519 1522 34.1 0.1532 0.1709 0.1529
> 1519 1521 32.5 0.1091 0.1221 0.1090
> 1519 1520 35.2 0.1091 0.1232 0.1090
> 1516 1519 38.5 0.1533 0.1791 0.1529
> 1516 1518 51.8 0.1093 0.1334 0.1090
> 1516 1517 56.4 0.1094 0.1351 0.1090
> 1513 1515 58.5 0.1092 0.1557 0.1090
> 1513 1514 56.9 0.1092 0.1569 0.1090
> 1510 1513 60.5 0.1545 0.1522 0.1529
> 1510 1512 53.6 0.1103 0.1148 0.1090
> 1510 1511 41.7 0.1101 0.1016 0.1090
> 1507 1510 55.0 0.1545 0.1515 0.1529
> 1507 1509 58.9 0.1093 0.1563 0.1090
> 1507 1508 54.7 0.1092 0.1557 0.1090
> 1504 1506 45.3 0.1093 0.1293 0.1090
> 1504 1505 49.1 0.1094 0.1310 0.1090
> 1501 1504 38.9 0.1533 0.1749 0.1529
> 1818 1820 34.4 0.1091 0.1088 0.1090
> 2563 2565 39.8 0.1093 0.1167 0.1090
> 2563 2564 37.6 0.1093 0.1163 0.1090
> 2560 2562 42.1 0.1092 0.1220 0.1090
> 2560 2561 40.0 0.1092 0.1214 0.1090
> 2557 2560 35.8 0.1536 0.1535 0.1529
> 2554 2557 33.7 0.1537 0.1525 0.1529
> 2554 2556 37.1 0.1092 0.1219 0.1090
> 2554 2555 40.0 0.1092 0.1226 0.1090
> 2551 2553 36.5 0.1093 0.1167 0.1090
> 2551 2552 52.9 0.1092 0.1149 0.1090
> Wrote pdb files with previous and current coordinates
> Step= 15, Dmax= 1.3e-01 nm, Epot= 4.73137e+04 Fmax= 1.84271e+06, atom=
> 822
> Step 16, time 0.032 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.005654, max 0.128942 (between atoms 825 and 828)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> 828 829 30.8 0.1113 0.1110 0.1090
> 825 828 31.7 0.1555 0.1726 0.1529
> 825 827 47.4 0.1102 0.1221 0.1090
> 825 826 76.4 0.1085 0.1166 0.1090
> 822 824 34.7 0.1092 0.1223 0.1090
> 1765 1767 37.2 0.1145 0.1126 0.1090
> 1765 1766 47.4 0.1128 0.1127 0.1090
> Step= 17, Dmax= 3.2e-02 nm, Epot= 3.09838e+03 Fmax= 2.83061e+03, atom=
> 1766
> Step= 18, Dmax= 3.9e-02 nm, Epot= 2.91458e+03 Fmax= 2.84186e+03, atom=
> 828
>
> Step 19, time 0.038 (ps) LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.003357, max 0.050861 (between atoms 825 and 826)
> bonds that rotated more than 30 degrees:
> atom 1 atom 2 angle previous, current, constraint length
> 828 830 39.1 0.1095 0.1118 0.1090
> 825 826 39.2 0.1087 0.1145 0.1090
> 1516 1517 31.3 0.1094 0.1133 0.1090
> Step= 20, Dmax= 2.3e-02 nm, Epot= 2.38983e+03 Fmax= 1.82748e+03, atom=
> 1516
> Step= 22, Dmax= 1.4e-02 nm, Epot= 2.05520e+03 Fmax= 1.32025e+03, atom=
> 1510
> Step= 23, Dmax= 1.7e-02 nm, Epot= 2.00275e+03 Fmax= 2.13164e+03, atom=
> 1516
> Step= 24, Dmax= 2.0e-02 nm, Epot= 1.86550e+03 Fmax= 2.34225e+03, atom=
> 1510
> Step= 25, Dmax= 2.4e-02 nm, Epot= 1.83379e+03 Fmax= 2.62066e+03, atom=
> 1519
> Step= 27, Dmax= 1.4e-02 nm, Epot= 1.41757e+03 Fmax= 5.65049e+02, atom=
> 1513
>
> writing lowest energy coordinates.
>
> Steepest Descents converged to Fmax < 1000 in 28 steps
> Potential Energy = 1.41757054835921e+03
> Maximum force = 5.65049499838719e+02 on atom 1513
> Norm of force = 7.42266925121288e+01
>
> ******************************
> *****************************************************************
>
> I tried a NVT for 500 ps and noticed PE LJ SR are reaching equilibrium
> (plateau after about 100ps). T coupling is working perfectly as well....
>
> Please let me know if my approach/results (picking a convoluted shape
> regardless of equilibrium state) are reliable...
>
> Thank you for your help :)
>
> Also I tried to compress a system of 8 polymers (rather than replicating a
> compressed single chain). As box size becomes less than 3.4 simulation
> crashes but what I need is 3.6nm so I got the desired density. A separate
> folder is included in attached file.(you may compare the results in 2 Excel
> files obtained from these two approaches)
>
> As far as my system is concerned I have no clear idea on the following
> issues:
>
> 1- In general in most of the NPT runs I am not reaching the reference
> pressure ( 50 ,30 ..bar). I get very low values around zero or 1.
> 2-About bond constraints: Actually, in literature I see it is common to fix
> the bond distance. Also in one of the GROMACS tutorials I have seen the
> reasoning that" since bond vibrations are of quantum nature they can not be
> treated classically". The reason you mentioned that I am noy using a good
> force filed is because of using bond constraints or ...? Do I need to remove
> constraints?
>
> 3- I have another inquiry about treating electrostatics contribution in my
> system. My systems as you have noticed involves only C and H. Later I will
> have to solvate this PE chain in HExane which is again apolar (only partial
> charges). It is not still clear to me whether I need to include
> electrostatic energy terms or I need to consider only vdw for calculating
> interaction energies
>
> 4-the value of compressibility I am using in mdp file is 4.5e-5 which is for
> water. I have not found this valued for my polymer yet.
>
> Thanks again, :)
>
>>
>> --
>
> Moeed Shahamat
> Graduate Student (Materials Modeling Research Group)
> McGill University- Department of Chemical Engineering
> Montreal, Quebec H3A 2B2, Canada
> Web:http://mmrg.chemeng.mcgill.ca/pages/current-group-members/moeed-shahamat.php
> Web:http://mmrg.chemeng.mcgill.ca/
>
More information about the gromacs.org_gmx-users
mailing list