[gmx-users] about my single point calculation
Hello Mark, I don't get any informing of your reply by e-mail, but get your reply searched by google. Anyway thanks very much for your reply! Yeah, I used two totally different mdp files for the single point calculation because I thought that gromcs would report the potential energy of my system if I used the option -rerun, no matter what mdp files I used. Gromacs-4.5.5 was used in the calculations. Some time later, I also tried as you mentioned in the e-mail below. Problem was the same. I used this mdp (attached minim.mdp) file for both 0-step minimization and single point calculation with rerun: The only difference is the integrator, steep for 0-step minimization while md for single point calculation. And the following lines are the output I got: Single point calculation with rerun : Step Time Lambda 00.00.0 Energies (kJ/mol) G96Bond G96Angle Improper Dih. LJ-14 Coulomb-14 3.97878e+051.44370e+041.06677e+043.93431e+01 9.36593e+01 LJ (SR)Coulomb (SR) Potential Kinetic En. Total Energy 2.77574e+015.17380e+024.23661e+050.0e+00 4.23661e+05 Temperature Pressure (bar) 0.0e+000.0e+00 0-step minimization: Step Time Lambda 00.00.0 Energies (kJ/mol) G96Bond G96AngleImproper Dih. LJ-14 Coulomb-14 5.75700e+011.78703e+018.25973e-02 -9.22001e+00 5.93500e+01 LJ (SR)Coulomb (SR) Potential Pressure (bar) -1.96671e+013.75963e+024.81949e+020.0e+00 Later, I also used the mdp file which was pasted on the forum (attached sp.mdp) to do single point calculation with rerun, The following is what I got: Step Time Lambda 00.00.0 Energies (kJ/mol) G96BondG96Angle Improper Dih. LJ-14 Coulomb-14 3.97878e+051.44370e+041.06677e+043.93431e+01 9.36593e+01 LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy 2.77575e+015.17379e+024.23661e+056.54617e+01 4.23726e+05 Conserved En.Temperature Pressure (bar) 4.23726e+051.45800e+020.0e+00 I found that those energies are pretty much the same as the one mentioned above. The differences between corresponding energies are huge, I still don't understand the difference. My system only contains 37 atom, the energies generated from the 0-step minimization seem more reasonable than those from single point calculation with rerun. Do you know the possible reasons which could result in the huge difference? Or I made some mistake? Thanks very much! All the best, Qinghua On Wed, Nov 6, 2013 at 4:07 PM, fantasticqhl fantastic...@gmail.com wrote: Dear Justin, I am sorry for the late reply. I still can't figure it out. It isn't rocket science - your two .mdp files describe totally different model physics. To compare things, change as few things as necessary to generate the comparison. So use the same input .mdp file for the MD vs EM single-point comparison, just changing the integrator line, and maybe unconstrained-start (I forget the details). And be aware of http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy Mark Could you please send me the mdp file which was used for your single point calculations. I want to do some comparison and then solve the problem. Thanks very much! All the best, Qinghua -- View this message in context: http://gromacs.5086.x6.nabble.com/single-point-calculation-with-gromacs-tp5012084p5012295.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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? Readhttp://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? Readhttp://www.gromacs.org/Support/Mailing_Lists define = ;-DPOSRES integrator = md ; molecular dynamics algorithm tinit = 0.0 ; start time and timestep in ps dt = 0.002; time step in ps nsteps = 2; number of steps for 1000ns run emtol = 100; convergence criterion emstep
[gmx-users] Re: single point calculation with gromacs
Dear Justin, I am sorry for the late reply. I still can't figure it out. Could you please send me the mdp file which was used for your single point calculations. I want to do some comparison and then solve the problem. Thanks very much! All the best, Qinghua -- View this message in context: http://gromacs.5086.x6.nabble.com/single-point-calculation-with-gromacs-tp5012084p5012295.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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: single point calculation with gromacs
Dear Justin, *Thanks very much for your reply! Here is my minim.mdp I used:* /; minim.mdp - used as input into grompp to generate em.tpr ; Parameters describing what to do, when to stop and what to save integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0; Stop minimization when the maximum force 1000.0 kJ/mol/nm emstep = 0.01 ; Energy step size nsteps = 0 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces ns_type = simple; Method to determine neighbor list (simple, grid) rlist = 9.0 ; Cut-off for making neighbor list (short range forces) coulombtype = Cut-off ; Treatment of long range electrostatic interactions rcoulomb= 9.0 ; Short-range electrostatic cut-off rvdw= 9.0 ; Short-range Van der Waals cut-off pbc = no; Periodic Boundary Conditions (yes/no) / *If I used this mdp file for 0-step minimization, I could get potential energy around 700 kJ/mol for my system. But I would get potential energy around 2.6e+05 kJ/mol if I used the mdv.mdp (md simulation in vacuum) for calculations with rerun option. And the bond potential could also reach as high as 1.1e+05 kJ/mol. Actually, my system is very small, it contains only 37 atoms. I believe that the energy reported by 0-step minimization were more reasonable. So I guess that there might be some problem for the mdp file usd with rerun, and here is the mdv.mdp:* /define = ;-DPOSRES integrator = md ; molecular dynamics algorithm tinit = 0.0 ; start time and timestep in ps dt = 0.002; time step in ps nsteps = 2; number of steps for 1000ns run emtol = 100; convergence criterion emstep = 0.05 ; intial step size nstlist = 0 ; step frequency for updating neighbour list ns_type = simple ; method for neighbour searching (?) nstxout = 1; frequency for writing coords to output .trr file nstvout = 1 ; frequency for writing velocities to output...should be same as nstxout nstfout = 1; frequency for writing forces to output nstlog = 1 ; frequency for writing energies to log file nstenergy = 1 ; frequency for writing energies to energy file nstxtcout = 1 ; frequency for writing coords to xtc traj xtc_grps= system ; group(s) whose coords are to be written in xtc traj energygrps = system ; group(s) whose energy is to be written in energy file pbc = no ; use pbc rlist = 0 ; cutoff lengths (nm) epsilon_r = 1.0 ; Dielectric constant (DC) for twin-range or DC of reaction field niter = 100 ; Some thingies for future use fourierspacing = 0.16 fourier_nx = 30 fourier_ny = 30 fourier_nz = 30 coulombtype = Cut-off ; truncation for minimisation, with large cutoff rcoulomb= 0 rcoulomb-switch = 0 vdw-type = Cut-off ; truncation for minimisation, with large cutoff rvdw-switch = 0 rvdw = 0 ; cut-off lengths epsilon_surface = 0 optimize_fft = yes Tcoupl = V-rescale tc_grps = system tau_t = 0.01 ref_t = 300 Pcoupl = no ; Parrinello-Rahman ; Pressure coupling gen_vel = yes gen_temp= 300 gen_seed= -1 constraints = none ; OPTIONS FOR BOND CONSTRAINTS constraint-algorithm = Lincs ; Type of constraint algorithm lincs_order = 4;4; Highest order in the expansion of the constraint coupling matrix lincs_iter = 1 lincs_warnangle = 30 ; Lincs will write a warning to the stderr if in one step a bond rotates ; over more degrees than unconstrained-start = no ; Do not constrain the start configuration ;Shake-SOR= no ; Use successive overrelaxation to reduce the number of shake iterations ;shake-tol= 1e-04 ; Relative tolerance of shake morse= no ; Convert harmonic bonds to morse potentials / *Could you please have a check for me again? Thanks in advance! All the best, Qinghua* -- View this message in context: http://gromacs.5086.x6.nabble.com/single-point-calculation-with-gromacs-tp5012084p5012100.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- gmx-users mailing list
[gmx-users] Re: single point calculation with gromacs
Dear Justin, *Thanks very much for your reply! Here is my minim.mdp I used:* /; minim.mdp - used as input into grompp to generate em.tpr ; Parameters describing what to do, when to stop and what to save integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0; Stop minimization when the maximum force 1000.0 kJ/mol/nm emstep = 0.01 ; Energy step size nsteps = 0 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces ns_type = simple; Method to determine neighbor list (simple, grid) rlist = 9.0 ; Cut-off for making neighbor list (short range forces) coulombtype = Cut-off ; Treatment of long range electrostatic interactions rcoulomb= 9.0 ; Short-range electrostatic cut-off rvdw= 9.0 ; Short-range Van der Waals cut-off pbc = no; Periodic Boundary Conditions (yes/no) / *If I used this mdp file for 0-step minimization, I could get potential energy around 700 kJ/mol for my system. But I would get potential energy around 2.6e+05 kJ/mol if I used the mdv.mdp (md simulation in vacuum) for calculations with rerun option. And the bond potential could also reach as high as 1.1e+05 kJ/mol. Actually, my system is very small, it contains only 37 atoms. I believe that the energy reported by 0-step minimization were more reasonable. So I guess that there might be some problem for the mdp file usd with rerun, and here is the mdv.mdp:* /define = ;-DPOSRES integrator = md ; molecular dynamics algorithm tinit = 0.0 ; start time and timestep in ps dt = 0.002; time step in ps nsteps = 2; number of steps for 1000ns run emtol = 100; convergence criterion emstep = 0.05 ; intial step size nstlist = 0 ; step frequency for updating neighbour list ns_type = simple ; method for neighbour searching (?) nstxout = 1; frequency for writing coords to output .trr file nstvout = 1 ; frequency for writing velocities to output...should be same as nstxout nstfout = 1; frequency for writing forces to output nstlog = 1 ; frequency for writing energies to log file nstenergy = 1 ; frequency for writing energies to energy file nstxtcout = 1 ; frequency for writing coords to xtc traj xtc_grps= system ; group(s) whose coords are to be written in xtc traj energygrps = system ; group(s) whose energy is to be written in energy file pbc = no ; use pbc rlist = 0 ; cutoff lengths (nm) epsilon_r = 1.0 ; Dielectric constant (DC) for twin-range or DC of reaction field niter = 100 ; Some thingies for future use fourierspacing = 0.16 fourier_nx = 30 fourier_ny = 30 fourier_nz = 30 coulombtype = Cut-off ; truncation for minimisation, with large cutoff rcoulomb= 0 rcoulomb-switch = 0 vdw-type = Cut-off ; truncation for minimisation, with large cutoff rvdw-switch = 0 rvdw = 0 ; cut-off lengths epsilon_surface = 0 optimize_fft = yes Tcoupl = V-rescale tc_grps = system tau_t = 0.01 ref_t = 300 Pcoupl = no ; Parrinello-Rahman ; Pressure coupling gen_vel = yes gen_temp= 300 gen_seed= -1 constraints = none ; OPTIONS FOR BOND CONSTRAINTS constraint-algorithm = Lincs ; Type of constraint algorithm lincs_order = 4;4; Highest order in the expansion of the constraint coupling matrix lincs_iter = 1 lincs_warnangle = 30 ; Lincs will write a warning to the stderr if in one step a bond rotates ; over more degrees than unconstrained-start = no ; Do not constrain the start configuration ;Shake-SOR= no ; Use successive overrelaxation to reduce the number of shake iterations ;shake-tol= 1e-04 ; Relative tolerance of shake morse= no ; Convert harmonic bonds to morse potentials / *Could you please have a check for me again? Thanks in advance! All the best, Qinghua* -- View this message in context: http://gromacs.5086.x6.nabble.com/single-point-calculation-with-gromacs-tp5012084p5012101.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- gmx-users mailing list
[gmx-users] single point calculation with gromacs
Dear All, I want to do the single point calculations for my systems with gromacs, I used the method mentioned on the gmx's website. mdrun -s sp.tpr -rerun configuration.pdb. However, I have some questions on how to generate the tpr file for single point calculation. It didn't work if I used a mdp files of minimization to generate the tpr file. It worked if I used a mdp files of md in vacuum, but the energies were much higher than those reported by the 0-step minimization. I guess that there might be some problem. Could some one give me some suggestion on the mdp file for this kind of single point callculations? Thanks very much! All the best, Qinghua -- View this message in context: http://gromacs.5086.x6.nabble.com/single-point-calculation-with-gromacs-tp5012084.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Energies in MM scanning using gmx
Dear all, But I have questions about the energies calculated by GMX for bonds and angles scanning. For example, a series of conformations of one system should be generated for a bond scanning, the only difference in geometry of these conformations is the bond length A---B. However, when I set the force constant as 0 for bond A---B, the bond energies is not a flat line as the bond length changes, in principle, it should be a flat line, but it is not. http://gromacs.5086.x6.nabble.com/file/n5010274/bond_1FF.png This is an example I got for my system. Could some tell me how to explain the energy differences? Thanks very much! All the best, qh -- View this message in context: http://gromacs.5086.x6.nabble.com/Energies-in-MM-scanning-using-gmx-tp5010274.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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: Energies in MM scanning using gmx
Thank you, all! I know the reason now! All the best, qh -- View this message in context: http://gromacs.5086.x6.nabble.com/Energies-in-MM-scanning-using-gmx-tp5010274p5010278.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Charges for calculation of Coulomb interaction
Dear GMX users, Just a simple question, are the Coulomb interactions calculated using charges from ffnonbonded.itp or aminoacid.rtp? Normally, they are the same. However, if I modify some charges for some atoms (for example CB) of some residues in the aminoacid.rtp, but do not modify in ffnonbonded.itp because some other residues would still used the old charges, how will gromacs treat the calculations? Right now I know the best option is introducing new atom types. Thanks very much! All the best, Qinghua -- 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: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks very much for your explanation. I will try to use the method as you suggested to do the validation. Thanks! All the best, Qinghua Liao On 04/09/2013 12:10 PM, Dr. Vitaly Chaban wrote: On Tue, Apr 9, 2013 at 11:03 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much again. I am sorry for the unclear, charge transfer was also taken into account for the complex, I did not mentioned in the last e-mail. What do you mean by finite T effect in MD? Kinetics? I mean thermal motion. You have an optimal structure/energy at 0K in QM. In MD you want to simulate at higher T, I guess. The optimal structures in both cases may be very similar, but may be not so similar. I am just saying that you should not expect ideal coincidence of energy vs. distance curves. For the reproduction of binding energy, I guess I know how to do it using QM method. Simply, I just need to do three single point calculations for complex, ligands and ion, respectively. And correct for BSSE. For MM method, it is similar, however, I am not sure I get get the MM energy for just one ion. This energy is zero within classical MD, since you do not consider electrons and nucleus, as you do in QM. Only one calculation is needed for MM. You define the charge groups, such as ion and ligand and look at the interaction between them (g_energy). Is my understanding right? Thanks for all your explanations and suggestions on this problems! All the best, Qinghua Liao On 04/09/2013 10:03 AM, Dr. Vitaly Chaban wrote: On Tue, Apr 9, 2013 at 9:39 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your patient and detailed suggestions on this problem. Actually, I am doing what your suggested now. I optimized the copper-ligand complex using QM method, and then did some QM scannings to derive the bond and angle force constants. Right now, I am doing the MM scanning using the same coordinates which were used in the QM scanning. What we want is that the MM curves can reproduce the QM curves. I think it is simply impossible in your case to reproduce the QM curves. You neglect charge transfer from copper to the ligand, resulting a chemical bond formation, you neglect finite T effect in your MD. If you want to remain in the framework of LJ+Coulomb, the best think you can get is reproduction of ion-ligand binding energy and more or less adequate distance ion-closest atom of the ligand But some of them agreed well, some of them did not. So I try to tune the sigma of the liganded atoms, however, it is a little complicated to tune many liganded atoms at the same time. I am still trying to work it out. Start from the sigma for ion-closest atom of the ligand. All other atoms will adjust automatically, since they are connected all together within the ligand. My personal viewpoint, which you may share or not, is not to do anything with sigmas of other atoms of the ligand. It is best for future portability to limit refinement to the ion only. It seems that you have much experience on such problems, could you please give me some suggestions on tuning the sigmas of atoms again? Thanks very much in advance! All the best, Qinghua Liao On 04/08/2013 03:51 PM, Dr. Vitaly Chaban wrote: On Mon, Apr 8, 2013 at 3:36 PM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your patient explanation. Yeah, you are right, that is what I want to know: how you tuned this parameter? Since then, if I want to set a new atom type and I know its vdw radius, so how should I set the sigma for it based on the vdw radius, You cannot set the sigma based ONLY on the VDW radius. which should be in agreement with OPLS-AA/L force filed? Could you give me some suggestions? I guess that I have to tune it by myself this time, right? Thanks in advance! I would do the following: 1) Optimize ion-ligand complex using ab initio. Write down binding energy and optimal distance; 2) Construct topology for classical MD using approximate sigma; 3) Calculate energy and distance from classical MD; 4) Compare them to distance and energy from ab initio; 5) If you are not satisfied, adjust your sigma; 6) Repeat classical MD until the difference between ion-ligand distance in classical MD becomes reasonably
Re: [gmx-users] Re: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks very much for your patient and detailed suggestions on this problem. Actually, I am doing what your suggested now. I optimized the copper-ligand complex using QM method, and then did some QM scannings to derive the bond and angle force constants. Right now, I am doing the MM scanning using the same coordinates which were used in the QM scanning. What we want is that the MM curves can reproduce the QM curves. But some of them agreed well, some of them did not. So I try to tune the sigma of the liganded atoms, however, it is a little complicated to tune many liganded atoms at the same time. I am still trying to work it out. It seems that you have much experience on such problems, could you please give me some suggestions on tuning the sigmas of atoms again? Thanks very much in advance! All the best, Qinghua Liao On 04/08/2013 03:51 PM, Dr. Vitaly Chaban wrote: On Mon, Apr 8, 2013 at 3:36 PM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your patient explanation. Yeah, you are right, that is what I want to know: how you tuned this parameter? Since then, if I want to set a new atom type and I know its vdw radius, so how should I set the sigma for it based on the vdw radius, You cannot set the sigma based ONLY on the VDW radius. which should be in agreement with OPLS-AA/L force filed? Could you give me some suggestions? I guess that I have to tune it by myself this time, right? Thanks in advance! I would do the following: 1) Optimize ion-ligand complex using ab initio. Write down binding energy and optimal distance; 2) Construct topology for classical MD using approximate sigma; 3) Calculate energy and distance from classical MD; 4) Compare them to distance and energy from ab initio; 5) If you are not satisfied, adjust your sigma; 6) Repeat classical MD until the difference between ion-ligand distance in classical MD becomes reasonably similar to that in ab initio. To preserve compatibility with OPLS, use the same level of theory in ab initio, which they used when derived OPLS. Keep in mind that their original level of theory is not so perfect... Dr. Vitaly Chaban All the best, Qinghua Liao -- 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: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks very much again. I am sorry for the unclear, charge transfer was also taken into account for the complex, I did not mentioned in the last e-mail. What do you mean by finite T effect in MD? Kinetics? For the reproduction of binding energy, I guess I know how to do it using QM method. Simply, I just need to do three single point calculations for complex, ligands and ion, respectively. For MM method, it is similar, however, I am not sure I get get the MM energy for just one ion. Is my understanding right? Thanks for all your explanations and suggestions on this problems! All the best, Qinghua Liao On 04/09/2013 10:03 AM, Dr. Vitaly Chaban wrote: On Tue, Apr 9, 2013 at 9:39 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your patient and detailed suggestions on this problem. Actually, I am doing what your suggested now. I optimized the copper-ligand complex using QM method, and then did some QM scannings to derive the bond and angle force constants. Right now, I am doing the MM scanning using the same coordinates which were used in the QM scanning. What we want is that the MM curves can reproduce the QM curves. I think it is simply impossible in your case to reproduce the QM curves. You neglect charge transfer from copper to the ligand, resulting a chemical bond formation, you neglect finite T effect in your MD. If you want to remain in the framework of LJ+Coulomb, the best think you can get is reproduction of ion-ligand binding energy and more or less adequate distance ion-closest atom of the ligand But some of them agreed well, some of them did not. So I try to tune the sigma of the liganded atoms, however, it is a little complicated to tune many liganded atoms at the same time. I am still trying to work it out. Start from the sigma for ion-closest atom of the ligand. All other atoms will adjust automatically, since they are connected all together within the ligand. My personal viewpoint, which you may share or not, is not to do anything with sigmas of other atoms of the ligand. It is best for future portability to limit refinement to the ion only. It seems that you have much experience on such problems, could you please give me some suggestions on tuning the sigmas of atoms again? Thanks very much in advance! All the best, Qinghua Liao On 04/08/2013 03:51 PM, Dr. Vitaly Chaban wrote: On Mon, Apr 8, 2013 at 3:36 PM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your patient explanation. Yeah, you are right, that is what I want to know: how you tuned this parameter? Since then, if I want to set a new atom type and I know its vdw radius, so how should I set the sigma for it based on the vdw radius, You cannot set the sigma based ONLY on the VDW radius. which should be in agreement with OPLS-AA/L force filed? Could you give me some suggestions? I guess that I have to tune it by myself this time, right? Thanks in advance! I would do the following: 1) Optimize ion-ligand complex using ab initio. Write down binding energy and optimal distance; 2) Construct topology for classical MD using approximate sigma; 3) Calculate energy and distance from classical MD; 4) Compare them to distance and energy from ab initio; 5) If you are not satisfied, adjust your sigma; 6) Repeat classical MD until the difference between ion-ligand distance in classical MD becomes reasonably similar to that in ab initio. To preserve compatibility with OPLS, use the same level of theory in ab initio, which they used when derived OPLS. Keep in mind that their original level of theory is not so perfect... Dr. Vitaly Chaban All the best, Qinghua Liao -- 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: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks very much for concern on my research! We are going to the use the bonded model together with Coulomb and LJ potentials. My problem is that vdw radius and its sigma do not follow the equation of Rvdw = pow(2, 1/6)*sigma in the OPLS force field files, not just for copper. That's why I sent these e-mails for suggestions. I am sorry for the unclear. All the best, Qinghua Liao On 04/08/2013 01:22 PM, Dr. Vitaly Chaban wrote: Dear Qinghua Liao - In that case, I am just wishing you luck with the copper containing systems. Are you going to simulate copper-ligand interactions using Coulomb+LJ potential only? I would guess it is a chemical bonding case. Maybe the Morse potential (additionally) can be of better service? Dr. Vitaly Chaban On Mon, Apr 8, 2013 at 1:09 PM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your explanation. I guess that I get what you mean now! Thanks! All the best, Qinghua Liao On 04/07/2013 11:35 AM, Dr. Vitaly Chaban wrote: The equation is a direct consequence of LJ-12-6 equation. This equation is used in OPLS and most other force fields. The difference you found originate from the fact that, besides LJ potential, there is much stronger Coulomb potential in the copper-ion case. If you run simulations, you will see that copper-ligand distance is smaller than the sum of their sigmas multiplied by pow (2, 1/6). Dr. Vitaly Chaban On Sun, Apr 7, 2013 at 11:28 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks for the explanation. I know this equation. However, the van der Waals radius and its counterpart sigma in OPLS-AA/L force field files do not follow this equation. For example, the vdw radius of copper ion is 1.4 angstrom, and its sigma is 2.08470e-01 (I guess the unit is nm). pow(2, 1/6) is more than 1, so obviously this equation does not work with copper. So do other atoms. I guess that there might be an additional coefficient for this equation in gromacs. That's the purpose for asking. Thanks very much! All the best, Qinghua On 04/07/2013 10:48 AM, Dr. Vitaly Chaban wrote: Dear Qinghua - The formal relation is diameter = pow (2, 1/6) * sigma, provided that you have only LJ potential in your interacting subsystem. If this is not the case, an optimal sigma can only be found iteratively. Dr. Vitaly Chaban On Sun, Apr 7, 2013 at 10:36 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your reply. My question is the relationship between van der Waals radius and sigma in the OPLS-AA/L force filed files of Gromacs. Of course I did ab initio optimizations of my system, but I do not know there is some relation between the optimal bond length (copper--atom of the ligand) and sigma. Could you please be more clear and give a little detailed explanation? Thanks very much! All the best, Qinghua On 04/06/2013 06:07 PM, Dr. Vitaly Chaban wrote: In systems of such kind, everything will depend on the atom of the ligand, which coordinated by copper ion. Perform ab initio geometry optimization and find the optimal distance. Then adjust sigma(s). Dr. Vitaly Chaban There is a copper ion with four ligands in my system. I am going to study this system using MD simulations. For the vdW parameters, R*=1.74 angstrom and epsilon=1.14 kcal.mol from one paper will be used in our simulations. I already found the parameters of copper ion (Cu2+) in the OPLS-AA/L force field files: sigma= 2.08470e-01, epsilon=4.76976e+00, which are for Cu2+ without ligands. The two epsilon are the same, just with different units. My question is that I do not know how to convert the vdW radius to sigma. I found that the vdw radius of copper is 1.4 angstrom, and the sigma in the force field file is 2.08470e-01. Could someone tell me how to do the converting? Thanks very much! -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org
Re: [gmx-users] Re: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks very much for your patient explanation. Yeah, you are right, that is what I want to know: how you tuned this parameter? Since then, if I want to set a new atom type and I know its vdw radius, so how should I set the sigma for it based on the vdw radius, which should be in agreement with OPLS-AA/L force filed? Could you give me some suggestions? I guess that I have to tune it by myself this time, right? Thanks in advance! All the best, Qinghua Liao On 04/08/2013 03:21 PM, Dr. Vitaly Chaban wrote: I think your misunderstanding comes from the belief that sigma (as they are tabulated in the force field files) should *exactly correspond* to the VDW diameter, as in encyclopedia. This is simply not the case. In reality, sigmas in the force fields are tuned in order to give right interatomic distances AFTER you turn on all the necessary potentials (Coulombic attraction in case of OPLS). Dr. Vitaly Chaban On Mon, Apr 8, 2013 at 3:14 PM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for concern on my research! We are going to the use the bonded model together with Coulomb and LJ potentials. My problem is that vdw radius and its sigma do not follow the equation of Rvdw = pow(2, 1/6)*sigma in the OPLS force field files, not just for copper. That's why I sent these e-mails for suggestions. I am sorry for the unclear. All the best, Qinghua Liao On 04/08/2013 01:22 PM, Dr. Vitaly Chaban wrote: Dear Qinghua Liao - In that case, I am just wishing you luck with the copper containing systems. Are you going to simulate copper-ligand interactions using Coulomb+LJ potential only? I would guess it is a chemical bonding case. Maybe the Morse potential (additionally) can be of better service? Dr. Vitaly Chaban On Mon, Apr 8, 2013 at 1:09 PM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your explanation. I guess that I get what you mean now! Thanks! All the best, Qinghua Liao On 04/07/2013 11:35 AM, Dr. Vitaly Chaban wrote: The equation is a direct consequence of LJ-12-6 equation. This equation is used in OPLS and most other force fields. The difference you found originate from the fact that, besides LJ potential, there is much stronger Coulomb potential in the copper-ion case. If you run simulations, you will see that copper-ligand distance is smaller than the sum of their sigmas multiplied by pow (2, 1/6). Dr. Vitaly Chaban On Sun, Apr 7, 2013 at 11:28 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks for the explanation. I know this equation. However, the van der Waals radius and its counterpart sigma in OPLS-AA/L force field files do not follow this equation. For example, the vdw radius of copper ion is 1.4 angstrom, and its sigma is 2.08470e-01 (I guess the unit is nm). pow(2, 1/6) is more than 1, so obviously this equation does not work with copper. So do other atoms. I guess that there might be an additional coefficient for this equation in gromacs. That's the purpose for asking. Thanks very much! All the best, Qinghua On 04/07/2013 10:48 AM, Dr. Vitaly Chaban wrote: Dear Qinghua - The formal relation is diameter = pow (2, 1/6) * sigma, provided that you have only LJ potential in your interacting subsystem. If this is not the case, an optimal sigma can only be found iteratively. Dr. Vitaly Chaban On Sun, Apr 7, 2013 at 10:36 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your reply. My question is the relationship between van der Waals radius and sigma in the OPLS-AA/L force filed files of Gromacs. Of course I did ab initio optimizations of my system, but I do not know there is some relation between the optimal bond length (copper--atom of the ligand) and sigma. Could you please be more clear and give a little detailed explanation? Thanks very much! All the best, Qinghua On 04/06/2013 06:07 PM, Dr. Vitaly Chaban wrote: In systems of such kind, everything will depend on the atom
Re: [gmx-users] Re: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks very much for your reply. My question is the relationship between van der Waals radius and sigma in the OPLS-AA/L force filed files of Gromacs. Of course I did ab initio optimizations of my system, but I do not know there is some relation between the optimal bond length (copper--atom of the ligand) and sigma. Could you please be more clear and give a little detailed explanation? Thanks very much! All the best, Qinghua On 04/06/2013 06:07 PM, Dr. Vitaly Chaban wrote: In systems of such kind, everything will depend on the atom of the ligand, which coordinated by copper ion. Perform ab initio geometry optimization and find the optimal distance. Then adjust sigma(s). Dr. Vitaly Chaban There is a copper ion with four ligands in my system. I am going to study this system using MD simulations. For the vdW parameters, R*=1.74 angstrom and epsilon=1.14 kcal.mol from one paper will be used in our simulations. I already found the parameters of copper ion (Cu2+) in the OPLS-AA/L force field files: sigma= 2.08470e-01, epsilon=4.76976e+00, which are for Cu2+ without ligands. The two epsilon are the same, just with different units. My question is that I do not know how to convert the vdW radius to sigma. I found that the vdw radius of copper is 1.4 angstrom, and the sigma in the force field file is 2.08470e-01. Could someone tell me how to do the converting? Thanks very much! -- 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: How to set the sigma and epsilon for Cu2+ in OPLS-AA/L
Dear Dr. Vitaly Chaban, Thanks for the explanation. I know this equation. However, the van der Waals radius and its counterpart sigma in OPLS-AA/L force field files do not follow this equation. For example, the vdw radius of copper ion is 1.4 angstrom, and its sigma is 2.08470e-01 (I guess the unit is nm). pow(2, 1/6) is more than 1, so obviously this equation does not work with copper. So do other atoms. I guess that there might be an additional coefficient for this equation in gromacs. That's the purpose for asking. Thanks very much! All the best, Qinghua On 04/07/2013 10:48 AM, Dr. Vitaly Chaban wrote: Dear Qinghua - The formal relation is diameter = pow (2, 1/6) * sigma, provided that you have only LJ potential in your interacting subsystem. If this is not the case, an optimal sigma can only be found iteratively. Dr. Vitaly Chaban On Sun, Apr 7, 2013 at 10:36 AM, fantasticqhl fantastic...@gmail.com mailto:fantastic...@gmail.com wrote: Dear Dr. Vitaly Chaban, Thanks very much for your reply. My question is the relationship between van der Waals radius and sigma in the OPLS-AA/L force filed files of Gromacs. Of course I did ab initio optimizations of my system, but I do not know there is some relation between the optimal bond length (copper--atom of the ligand) and sigma. Could you please be more clear and give a little detailed explanation? Thanks very much! All the best, Qinghua On 04/06/2013 06:07 PM, Dr. Vitaly Chaban wrote: In systems of such kind, everything will depend on the atom of the ligand, which coordinated by copper ion. Perform ab initio geometry optimization and find the optimal distance. Then adjust sigma(s). Dr. Vitaly Chaban There is a copper ion with four ligands in my system. I am going to study this system using MD simulations. For the vdW parameters, R*=1.74 angstrom and epsilon=1.14 kcal.mol from one paper will be used in our simulations. I already found the parameters of copper ion (Cu2+) in the OPLS-AA/L force field files: sigma= 2.08470e-01, epsilon=4.76976e+00, which are for Cu2+ without ligands. The two epsilon are the same, just with different units. My question is that I do not know how to convert the vdW radius to sigma. I found that the vdw radius of copper is 1.4 angstrom, and the sigma in the force field file is 2.08470e-01. Could someone tell me how to do the converting? Thanks very much! -- 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 set the sigma and epsilon for Cu2+ in OPLS-AA/L force field
Dear GMX users, There is a copper ion with four ligands in my system. I am going to study this system using MD simulations. For the vdW parameters, R*=1.74 angstrom and epsilon=1.14 kcal.mol from one paper will be used in our simulations. I already found the parameters of copper ion (Cu2+) in the OPLS-AA/L force field files: sigma= 2.08470e-01, epsilon=4.76976e+00, which are for Cu2+ without ligands. The two epsilon are the same, just with different units. My question is that I do not know how to convert the vdW radius to sigma. I found that the vdw radius of copper is 1.4 angstrom, and the sigma in the force field file is 2.08470e-01. Could someone tell me how to do the converting? Thanks very much! All the best, Qinghua -- 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 modify the vdw radii of some atom?
Hi all, I want to modify the vdw radii of some atoms, but I did not find the file that I can edit for this purpose. I checked the file vdwradii.dat, but I did not find the atoms that I want to modify. Could someone tell me where I can modify them? Thanks very much! All the best, -- View this message in context: http://gromacs.5086.n6.nabble.com/How-to-modify-the-vdw-radii-of-some-atom-tp5006497.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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
