Hi. I’m using the LAMMPS colvars module to evaluate the free energy of adsorption of a molecule from a solid surface. I have run some simple tests from the example script in github and applied the module to my case, in which I found some questions. For all the questions, I have attached the colvars configuration file and the simulated videos in this link: https://drive.google.com/drive/folders/1MnsqPZKNrSivFZLo2UxvVmrO2bAL8KO5?usp=sharing . I will include the colvars content here as well.
- Question 1: For this testing case, a water molecule pulled away from the solid calcite surface in Z direction, I used the distanceZ collective variable and umbrella sampling module, setting the center 5 to target center 15. The molecule is pulled awayed in the intended direction, Z axis but the molecule suddenly fluctuates around in the Z direction… May I know what I should implement in the code to ensure it is pulled consistently without moving around? The PMF result generated looks logical anyway.
colvar {
name dist
width 0.1 # Spacing of grids
lowerBoundary 0.0
upperBoundary 15.0 # Determined by size of unit cell used
hardLowerBoundary yes # Inform biases that lower boundary won’t be crossed
hardUpperBoundary yes # Inform biases that upper boundary won’t be crossed
outputAppliedForce on # keep track of bias force on this variable
distanceZ {
axis (0.0, 0.0, 1.0) #Ensure the projection of distance vector on z axis
main {
atomNumbersRange 301-303
}
ref {
atomNumbersRange 1-300
}
}
}
harmonic {
name h_pot
colvars dist
forceConstant 1
centers 5
targetCenters 15
targetNumStages 19
targetNumSteps 50000
outputCenters yes # Write the current centers to the trajectory file
writeTIPMF yes
#outputAccumulatedWork yes # keep track of how much energy will be added
}
- Question 2: For the adaptive biasing force abf module. The abf module looks easier to run with the simple code. When I visualize the file, the molecule is just maintained in the same initial position till the end of the simulation. However, the PMF module shows the grid spacing from the initial distance to the final distance, and the PMF value is the same along the coordinate. The results seem weird and how should I use ABF properly?
colvar {
name dist
width 0.1 # Spacing of grids
lowerBoundary 0.0
upperBoundary 15.0 # Determined by size of unit cell used
hardLowerBoundary yes # Inform biases that lower boundary won’t be crossed
hardUpperBoundary yes # Inform biases that upper boundary won’t be crossed
outputAppliedForce on # keep track of bias force on this variable
distance {
group1 {
atomNumbersRange 1-300
}
group2 {
atomNumbersRange 301-303
}
}
}
abf {
colvars dist
fullsamples 200
maxForce 100
}
- Questions 3: I have also tried to modify the example script for umbrella sampling method - peptide run from this link in github: https://github.com/CFDEMproject/LAMMPS/tree/master/examples/USER/colvars. The simulation can run, but the reference atom 37 keeps moving together with the moving atom even though I have applied centre restraint to the atom with respect to the dummy atom. I have also tried another method from github using distanceVec and apply the harmonic restraint between dummy atom and reference atom, however the reference atom 37 will still move together… How can I hold the reference atom in the fixed position?
colvar {
name d
width 0.1 # typical displacement is set to 1 Angstrom
the position restraint is enforced as a distance to a fixed point
defined by a dummy atom group.
distance {
group1 {
group definition:
atomNumbers 37
centerReference # use relative coordinates
refPositionsFile data.peptide # initial coordinates for reference group # atom 37 position 43.51012 53.02289 44.10510
refPositions (43.51012, 53.02289, 44.10510)
}
group2 {
dummyAtom (0.0, 0.0, 0.0) # fixed position
}
}
}
harmonic { # Define a moving harmonic restraint
colvars d # acting on d colvars
centers 0 # centered around zero displacement, zero RMSD, and no rotation
forceConstant 10 # unit is kcal/mol/[width]^2
where the width parameter is adapted to the dimension of each colvar
}
colvar {
name dist
width 0.1 # Spacing of grids
lowerBoundary 0.0
upperBoundary 12.0 # Determined by size of unit cell used
hardLowerBoundary yes # Inform biases that lower boundary won’t be crossed
hardUpperBoundary yes # Inform biases that upper boundary won’t be crossed
outputAppliedForce on # keep track of bias force on this variable
distance {
forceNoPBC yes
group1 {
atomNumbers 2 4 5 6
}
group2 {
atomNumbers 37
}
}
}
harmonic {
name h_pot
colvars dist
forceConstant 1.0
centers 4.0
targetCenters 12.0
targetNumStages 9
targetNumSteps 10000
outputCenters yes # Write the current centers to the trajectory file
writeTIPMF yes
outputAccumulatedWork yes # keep track of how much energy will be added
}
Thank you!