[gmx-users] constraint/restrain distance
ÿffffc1ngel Piÿfffff1eiro
anxomx at yahoo.es
Fri Aug 20 14:39:52 CEST 2004
Dear gromacs users,
I want to calculate the force to fix the distance
between two molecules A and B. I think that it is not
possible at the moment to restraint the distance
between different molecules in the topology, so I
tried to use the pull code. I started to play with the
3 different options (AFM, CF and umbrella sampling)
but I found several problems.
1- I have been trying to constraint the distance using
a pull.ppa file like this:
runtype = constraint
group_1 = B
reference_group = A
reftype = com
Skip steps = 10
pulldim = Y Y Y
constraint_direction = 0.0 0.0 0.0
constraint_rate = 0
constraint_tolerance = 1E-6
I obtained very noisy plots although the averages at
different distances seem to be reasonable (when B is
too close to A the force is very positive and when the
distance is slightly higher than the expected one the
force is negative). The problem was that the
simulations crashed after several nanoseconds, or
picoseconds in some cases. The last lines of the log
file are:
***************************************************************************************************************
103 106 86.7 1.2344 21318.3443 0.1520
104 105 95.4 0.4626 892.0601 0.1000
106 107 92.5 0.4006 21093.5763 0.1435
107 108 93.3 0.0650 3351.6941 0.1000
Constraint error in algorithm Lincs at step 1704010
t = 3408.020 ps: Water molecule starting at atom 7540
can not be settled.
Check for bad contacts and/or reduce the
timestep.Wrote pdb files with previous and current
coordinates
Grid: -2147483648 x -2147483648 x -2147483648 cells
Fatal error: ci = -1 should be in 0 .. -1 [FILE
nsgrid.c, LINE 218]
***************************************************************************************************************
So I think that this error has something to do with
Lincs...
I tried to avoid this error in different ways: running
the simulation in double precision, decreasing the
time step to 1 and 0.5 fs (the original timestep was 2
fs) and I also used different values of
constraint_tolerance. Decreasing the time step
produced less crashes and the force did not change
significatively at a given distance with any of these
parameters. My system is quite small since A and B
have 140 atoms in total (I need around 2500 water
molecules to solvate them) so it does not need much
time to equilibrate. Averaging through the time before
crashing, the forces seem to be very reasonable
(although the error is very big due to the noise) but
I am still worried about this crashings... Can anyone
help me?
(Both topologies have been created by hand and they
have been working fine for a lot of simulations, even
simulations at 1000 kelvins did not crash so I do not
think that it is a problem in the topologies...)
2.- I have also tried the umbrella sampling (in order
to compare the results obtained from the CF). My
pull.ppa file looks like:
verbose = yes
runtype = umbrella
group_1 = B
reference_group = A
reftype = com
reflag = 1
Skip steps = 10
pulldim = Y N N
k1 = 1000
pos1 = -0.3 0 0
I used only one coordinate because I just want to
restrain the relative distance but it did not work
(the distance between A and B changes a lot). Should I
use "pulldim = Y Y Y" ? In such a case the output will
be the time and the 3 coordinates (the deviation of B
from the restrained position (x,y,z) ). How are the
forces calculated in this case? I do not know if the
movement of each dimension is restrained
independently... Is F=-kd (with d=sqrt(x*x+y*y+z*z))?
3.- I am now trying AFM with afm_rate1=0 but I again
have problems with the dimensions. In the simulations
I am running now I am using the following pull.ppa:
verbose = yes
runtype = afm
group_1 = B
reference_group = A
reftype = com
reflag = 1
Skip steps = 1
pulldim = Y Y Y
afm_rate1 = 0
afm_k1 = 2000
afm_dir1 = 1 0 0
afm_init1 = 0.3 0 0
I obtained a pull.pdo file with 10 columns (the time
and the coordinates of A, B and the spring
respectively). I guess that the force can be
calculated just from the distance between the spring
and B and I created a simple script to calculate it
but the distances obtained are too big (around 3 nm).
This is the script:
awk '{if (NR > 7) print $1 " " $8-$5 " " $9-$6 " "
$10-$7}' ../pull.pdo |
awk '{print $1 " " sqrt($2*$2+$3*$3+$4*$4)}' >
dist.out
The initial configuration was created using editconf
and the position of B relative to A is (0.3, 0, 0) for
this particular example, so by choosing afm_init1 =
0.3 0 0 the coordinates of the spring and the
coordinates of B should be equal (I tried also with
afm_init1= -0.3 0 0) and the result was similar.
Moreover I calculated the distances between A and B
with g_dist and they are always < 0.5 nm.
Other thing that I do not understand is that the
coordinate y of the molecule A always is the same as
the coordinate x of the spring, and the coordinate z
of both, A and the spring, are also always the same.
Take a look to some lines of pull.pdo:
***************************************************************************************************************************************************************
161.630005 3.723072 4.007028
1.630497 3.366747 3.423072
3.723034 4.007028 1.724641
1.630497
161.632004 3.723068 4.007199
1.630563 3.366929 3.423068
3.723042 4.007199 1.724370
1.630563
161.634003 3.723071 4.007378
1.630636 3.367114 3.423071
3.723065 4.007378 1.724109
1.630636
161.636002 3.723074 4.007556
1.630714 3.367295 3.423074
3.723097 4.007556 1.723856
1.630714
161.638000 3.723071 4.007728
1.630795 3.367466 3.423071
3.723132 4.007728 1.723608
1.630795
161.640015 3.723063 4.007890
1.630878 3.367629 3.423063
3.723170 4.007890 1.723364
1.630878
161.642014 3.723051 4.008047
1.630965 3.367785 3.423051
3.723214 4.008047 1.723125
1.630965
161.644012 3.723038 4.008197
1.631058 3.367936 3.423038
3.723265 4.008197 1.722891
1.631058
161.646011 3.723024 4.008343
1.631156 3.368084 3.423024
3.723324 4.008343 1.722665
1.631156
***************************************************************************************************************************************************************
Sorry for so long message but I hope someone can help
me at least in some of these problems and that this
can be useful for other people...
I should also say that I have been using gmx3.2.1 for
all these calculations.
Angel Pineiro.
PS: I am not expert in gromacs and this is the first
time I am doing this kind of calculations... (just in
case somethings sound strange)
______________________________________________
Renovamos el Correo Yahoo!: ¡100 MB GRATIS!
Nuevos servicios, más seguridad
http://correo.yahoo.es
More information about the gromacs.org_gmx-users
mailing list