Hi everyone,
I’ve recently faced a problem with the statistical result of my simulation in which I mixed a single Acetonitrile (ACN) or CO2 molecule in a bulk of water ( N_w 2880 ), and I aimed to measure the solvation-free energy of the mixture at a fixed scaled Lj and Coloumb potential.
I’ve applied "fix adapt/fep " to adapt the lambda value and atoms’ charge of CO2 or ACN. I then used the Fep compute command with two options “volume yes tail yes” since the simulation is performed under NPT ensemble with the pressure 1 atm and temperature 300 K with the tail correction of LJ potential.
I’ve used Trappe force fields for ACN (united atom version) and CO2 and SPC/E and TIP4P/epsilon force fields for water.
I’ve got the wrong pressure for the intermediate lambda values, which is always more than the desired value of 1 atm.
To overcome this issue, I started with the LAMMPS example on LAMMPS/example/PACKAGES/FEP/CH4hyd/fep01, and I submitted this example for 2 M MD steps to evaluate both simulation details.
It then turned to me that the LAMMPS example also presents the wrong pressure.
I’ve attached one of my LAMMPS scripts as an instance for you, to follow up on what I meant by having a look at the script, and I would be grateful if you have any suggestions, clues or queries to solve this issue.
In addition, I’ve uploaded all other required files to run this simulation on your system if you wish.
Best regards,
Russell.
in.fepti.acn.h2o (4.3 KB)
C2H3N_UA.txt (488 Bytes)
init_conf.data (1.4 MB)
PARMS.lmp (2.0 KB)
This seems to be a very specific issue with the code in the FEP package. Too specific for the kind of issue that people responding here can comment on. Your best shot is to contact the (main) author of the FEP package. He is generally very responsive to questions from users of his package and his contact details should be the README file in the folder of the package in the LAMMPS sources.
1 Like
Great!
Thank you so much Alex!
Regarding the above issue
I’ve checked my script, removed the free energy computation part, and attempted to simulate a system of water molecules and a single ACN or CO2 solute molecule under the NPT ensemble. However, I kept both types of molecules as rigid bodies via two different algorithms, the Shake algorithm for the water molecules and the other molecules of ACN or CO2 are treated via either the “rigid/small” or the “rigid/nvt/small” fixes.
I’ve considered the shake command the last fix command before the “fix npt” as the LAMMPS guide emphasized, and I excluded the ACN molecule from the “fix npt” group-ID to avoid the double time integration process. This structure works just for \lambda either 0 or 1, and it presents an average pressure more than the desired pressure value ( ~ 1 atm ) for the intermediate states, which means \lambda \in (0,1).
I put my foot a bit further and examined various situations of the solvent and solute, and I realized the system exhibits consistent pressure on 1 atm (the expected value) if both or one of the system’s members is flexible.
So, I’ve come here with a question, why can not we combine two separated rigid algorithms in a system? Please be aware that I’ve applied the Shake method to water molecules and the Rigid method to ACN or CO2.
I would be grateful if you have any suggestions, clues or queries to solve this issue.
Best
Russell
in.npt.h2o.acn (3.9 KB)
Please note, that there is no benefit from using fix shake over fix rigid for the water molecules, when you need to use fix rigid for other molecules in the system. Those will require using a shorter timestep for non-divergent integration of the rotational equations of motion and thus the benefit of fix shake to allow a larger timestep is moot.
You input has:
pair_style hybrid &
lj/cut/coul/long 14.0 10.0 &
lj/cut/tip4p/long/soft 1 2 1 1 0.105 1 0.5 10.0 14.0 10.0 #############n =1 here
When you are using a tip4p style for coulomb in your system, you must use the same pair style (and not as a hybrid style) for all coulomb, otherwise the interactions between the TIP4P water and other coulomb is going to be incorrect as other pair styles will not determine and use the M point as charged site but use the oxygen location instead.
timestep 1.0
This timestep is too large for stable time integration for small molecules with fix rigid. You’ll need something like 0.25
You can, but you will need to use the suitable settings.