Hi everyone,
Lammps version: 3 June 2022
OS: Cent OS
I am trying to do some free energy calculations using FEP method using compute FEP command paired with the lj/cut/soft pair style in lammps. For this, I am running a simulation during which the \lambda parameter is decreased gradually after a fixed number of steps. And for each of the \lambda run, compute FEP command calculates the perturbation energy. This is accomplished using the following lines in input file:
variable simlength equal 20000
variable nsteps equal 5
variable nruns equal ${nsteps}+1
variable totalsimlength equal ${simlength}*${nruns}
variable samplefreq equal 10
variable l file Lambdas.txt
variable lambda equal next(l)
variable delta_lam_fwd equal 0.05
variable delta_lam_bwd equal -0.05
fix TVSTAT Rest nvt temp ${TK} ${TK} 100
fix FEP_vdW all adapt/fep ${simlength} pair lj/cut/soft lambda 1 1*6 v_lambda pair lj/cut/soft lambda 2 2*6 v_lambda scale no after yes
compute Fwd_FEP all fep ${TK} pair lj/cut/soft lambda 1 1*6 v_delta_lam_fwd pair lj/cut/soft lambda 2 2*6 v_delta_lam_fwd tail yes volume no
compute Bwd_FEP all fep ${TK} pair lj/cut/soft lambda 1 1*6 v_delta_lam_bwd pair lj/cut/soft lambda 2 2*6 v_delta_lam_bwd tail yes volume no
fix deltaU all ave/time 1 1 ${samplefreq} c_Fwd_FEP[1] c_Bwd_FEP[1] mode scalar file deltaUs_all.dat
thermo 1000
thermo_style custom step etotal ke pe evdwl temp vol press v_delta_lam_fwd v_delta_lam_bwd
So, what happens is, the simulation runs as expected till the second last \lambda value, but for the last one, the compute Fwd_FEP command starts to give values same as those for compute Bwd_FEP, which is caused by the value of the variable v_delta_lam_fwd changing to -0.05 instead of 0.05 as set in the input script, as printed using the thermo command:
Step TotEng KinEng PotEng E_vdwl Temp Volume Press v_delta_lam_fwd v_delta_lam_bwd
...
350093000 -86136.311 4273.3368 -90409.648 8735.9442 299.85692 75834.308 730.42089 0.05 -0.05
350094000 -86106.398 4284.4323 -90390.83 8918.2387 300.63548 75834.308 1234.6366 0.05 -0.05
350095000 -86135.691 4246.2832 -90381.974 8694.3414 297.95859 75834.308 510.14636 0.05 -0.05
350096000 -85976.684 4261.9146 -90238.598 8660.8715 299.05544 75834.308 496.38937 0.05 -0.05
350097000 -85919.287 4312.6694 -90231.956 8917.6138 302.61686 75834.308 1409.9869 0.05 -0.05
350098000 -86058.619 4247.5722 -90306.191 8840.1815 298.04904 75834.308 1079.8458 0.05 -0.05
350099000 -85896.704 4267.1954 -90163.899 8711.1602 299.42598 75834.308 820.38633 0.05 -0.05
350100000 -85985.651 4268.8908 -90254.542 8670.4515 299.54495 75834.308 545.48806 0.05 -0.05 # second last lambda run ends here.
350101000 -85966.823 4303.1725 -90269.995 8707.6271 301.95047 75834.308 706.52829 -0.05 -0.05
350102000 -85982.031 4254.3892 -90236.42 8755.4484 298.52738 75834.308 811.1872 -0.05 -0.05
350103000 -86021.924 4221.8595 -90243.783 8909.3352 296.24479 75834.308 1328.6266 -0.05 -0.05
350104000 -86018.547 4242.1872 -90260.734 8780.0258 297.67118 75834.308 821.86939 -0.05 -0.05
350105000 -86006.549 4258.5056 -90265.055 8685.4644 298.81623 75834.308 629.8689 -0.05 -0.05
350106000 -85897.769 4280.0067 -90177.775 8701.0648 300.32494 75834.308 635.17761 -0.05 -0.05
350107000 -86031.247 4208.3746 -90239.621 8684.3429 295.29857 75834.308 606.78614 -0.05 -0.05
350108000 -85894.196 4251.2076 -90145.404 8779.1583 298.30413 75834.308 906.35704 -0.05 -0.05
350109000 -85948.198 4193.3 -90141.498 8665.2021 294.2408 75834.308 456.46995 -0.05 -0.05
350110000 -85846.401 4235.3346 -90081.735 8811.3486 297.19033 75834.308 1267.8322 -0.05 -0.05
350111000 -85995.452 4256.681 -90252.133 8794.8207 298.68819 75834.308 946.96043 -0.05 -0.05
350112000 -85961.304 4272.6076 -90233.911 8781.2769 299.80575 75834.308 928.04201 -0.05 -0.05
350113000 -86201.837 4126.5484 -90328.385 8779.8829 289.55688 75834.308 835.88337 -0.05 -0.05
350114000 -85968.249 4196.3182 -90164.567 8655.3571 294.45258 75834.308 532.46323 -0.05 -0.05
350115000 -85897.515 4272.6101 -90170.125 8612.3793 299.80593 75834.308 487.07478 -0.05 -0.05
350116000 -85855.7 4341.0402 -90196.74 8691.2443 304.60762 75834.308 571.6355 -0.05 -0.05
350117000 -85963.257 4210.5771 -90173.834 8798.0102 295.45312 75834.308 999.98062 -0.05 -0.05
350118000 -86018.397 4240.9741 -90259.371 8609.4535 297.58606 75834.308 320.38252 -0.05 -0.05
350119000 -85891.813 4312.3972 -90204.21 8649.7025 302.59776 75834.308 633.75554 -0.05 -0.05
350120000 -85944.851 4302.84 -90247.691 8664.7604 301.92714 75834.308 608.32034 -0.05 -0.05
The same behaviour is seen even if I run the simulation for just 1 \lambda value. Can someone help me identify what might be causing this? And what mistake I may be making.
I am attaching the input files as well.
data.n60_nocharge_eqm.lmp (1.2 MB) # input restart file
deltaUs_all.dat (315.8 KB) # fep data file (see last few lines)
in.lmp (2.1 KB) # input script
log.lammps (28.3 KB) # simulation log file
pair.lmp (3.1 KB) # pair coefficient definitions
Lambdas.txt (65 Bytes) # File from which \lambda values are read using next command.
Thanks, and regards,
Chaitanya