Hi all,
I'm trying to modify the COM_DISTANCE colvar to only act along a specified
axis. I've been using Feng Pan's excellent tutorial (here
<
https://ambermd.org/tutorials/advanced/tutorial31/custom_cv/index.html> or
in nfe_colvar.F90) but I'm bumping up against the limits of my knowledge of
fortran and the code.
I decided to follow Feng and use cv_r to define which component(s) to use
for the bias. However, it seemed like the best way to implement the changes
would be to element-wise multiply the center of mass coordinates by the
vector provided in cv_r. I copied all of the COM_DISTANCE
functions/subroutines and added code to pull in and multiply the cv_r
vector by the COMs. I added "_AXIS" to the ends of the names to distinguish
them from the originals. The code compiles and runs fine and the values I
get align with the position of the biases. (i.e. SMD with the potential
moving from pos_1 to pos_2 appears to work)
The issue I'm having is that the CV values aren't accurate. As a control, I
did a 2D PMD simulation with COM_DISTANCE for one CV and my modified CV for
the second axis. All settings were otherwise identical. (and bias was 0 for
both) For the modified CV, I set cv_r to 1, 1, 1 which should use all three
axes and produce identical results to COM_DISTANCE. But it doesn't - the
values are very different. Furthermore, if I change cv_r to 0, 0, 1 (should
use only the z-coordinate), I get identical results to using 1, 1, 1. So
something in my logic and/or program is wrong but I can't find the error. I
put the first few lines of output using all three CVs at the end of the
message for illustration. (the 2nd column is COM_DISTANCE and the 3rd and
4th are my CV with 1, 1, 1 and 0, 0, 1; respectively) I also pulled the
coordinates out of the trajectory and calculated the distances manually and
COM_DISTANCE is correct while both instances of my CV are wrong.
I'm attaching my modified nfe_colvars here along with the mdin and colvars
file I'm using. I'm using prmtop and inpcrd.4 from $AMBERHOME/test/nfe as
inputs.
Thanks for the help!
Gard
# = NFE%PMD
==================================================================
# << anchor(1) : position = 5.000000, 15.000000, 15.000000, 25.000000
# strength = 0.000000, 0.000000 >>
# << anchor(2) : position = 5.000000, 15.000000, 15.000000, 25.000000
# strength = 0.000000, 0.000000 >>
# << anchor(3) : position = 5.000000, 15.000000, 15.000000, 25.000000
# strength = 0.000000, 0.000000 >>
#
-----------------------------------------------------------------------------
# MD time (ps), CV(1:3)
#
=============================================================================
0.0000 17.45140271 5.65492630 5.65492630
0.4000 17.87154781 5.78842885 5.78842885
0.8000 17.95033220 5.80990116 5.80990116
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
- application/octet-stream attachment: mdin
- application/octet-stream attachment: cv_pmd.in
Received on Fri Aug 16 2024 - 16:00:01 PDT
