[lammps-users] what does R0 mean in "fix smd"?

LAMMPS LAMMPS Mailing List Mirror
1

Hi Axel,

in general, please always ask these kind of questions
on the lammps list, so that people running into the
same problem can find the solution in the mailing list
archives. i have very little time and don't want to
answer the same questions over and over again. thanks.

I am still confused about the R0, I've write 2 version of code: fix

just forget about R0. it has nothing to do with R.
set it to zero and be happy.

smd cvel 1000 -0.001 tether NULL NULL 37 0 and fix smd cvel 1000
-0.0010 NULL NULL 0 37, what is the difference between them?

consider that you have a spring and a (strong) wire. the length
of the wire is R0. your group is connected to this wire and the
spring which is then connected to your reference point. the wire
always keeps the moving part of the spring at a distance of R0
from your group. this makes sense for fix spring, but not in
steered MD.

And here is part of doc:
In cvel mode the distance R is incremented or decremented monotonously
according to the pulling (or pushing) velocity.

the distance R is determined by the COM of a specified group and the
reference tether point. How can you increase it?

the sentence you quoted above answers your question.

And I've found a error in my in.script: no integrator is applied to
the sio2 group. Is it why the group is not moving?

yes.

2

Hi, Axel

I am still confused about what happened in fix smd. You said:

consider that you have a spring and a (strong) wire. the length
of the wire is R0. your group is connected to this wire and the
spring which is then connected to your reference point. the wire
always keeps the moving part of the spring at a distance of R0
from your group. this makes sense for fix spring, but not in
steered MD.

When R0 is set to be zero, the model is simplified like this:
sio2---spring---ref.point
  >______________|
               >
              R
Now I use constant velocity mode, set velocity to be 1.0. Which end of
the spring is moving, sio2 end or ref.point end?

#logout:
LAMMPS (18 Apr 2010)
# Created by charmm2lammps v1.8.1 on Sun Apr 4 20:05:41 CST 2010
boundary p p p
units real
neigh_modify delay 2 every 1

atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic

pair_style lj/charmm/coul/charmm 8 10
pair_modify mix arithmetic
#kspace_style pppm 1e-4

read_data data.test
Scanning data file ...
  4 = max bonds/atom
  6 = max angles/atom
  18 = max dihedrals/atom
  1 = max impropers/atom
Reading data file ...
  orthogonal box = (-24 -24 -33.6931) to (24 24 83.5)
  1 by 1 by 1 processor grid
  27324 atoms
  27324 velocities
  21428 bonds
  24562 angles
  31944 dihedrals
  144 impropers
Finding 1-2 1-3 1-4 neighbors ...
  4 = max # of 1-2 neighbors
  12 = max # of 1-3 neighbors
  16 = max # of 1-4 neighbors
  20 = max # of special neighbors

special_bonds charmm

#define groups
group sio2 type 17 18 19 20 21
354 atoms in group sio2
group water type 2 13
17610 atoms in group water
group other subtract all sio2
26970 atoms in group other
group lipid subtract other water
9360 atoms in group lipid
region lower block -INF INF -INF INF -INF -15 units box
group lower region lower
4526 atoms in group lower
group lower_lipid intersect lower lipid
1201 atoms in group lower_lipid

#minimize
#minimize 1.0e-10 1.0e-12 100 1000

# fix nvt on water and lipid
fix 1 other nvt temp 323.15 323.15 100
# fix all tip3p water model
fix 2 water shake 1e-6 500 0 b 17 a 31
Finding SHAKE clusters ...
  0 = # of size 2 clusters
  0 = # of size 3 clusters
  0 = # of size 4 clusters
  5870 = # of frozen angles

fix rigid sio2 rigid molecule
1 rigid bodies with 354 atoms
# hold lower lipid bilayer
fix hold_lipid lower_lipid spring/self 10
# pull
fix pull sio2 smd cvel 50 -0.001 tether NULL NULL 37 0
# fix 2_sio2 sio2 nve
fix sio2 sio2 setforce 0 0 NULL
# compute sio2 sio2 reduce sum fz
# output forces acted upon sio2 sphere
#fix f_sio2 sio2 ave/time 1 500 1000 f_pull[3] f_pull[4]
f_pull[5] f_pull[6] f_pull[7] file f_sio2.addH.xx
# fix all
# generate initial velocity
velocity other create 323.15 1468 dist uniform

velocity sio2 set 0 0 -0.001 units box
thermo 10
thermo_style custom step cpu temp press f_pull[3] f_pull[4]
f_pull[5] f_pull[6] f_pull[7]

# 1fs per step
timestep 1

dump 1 all atom 10 dump.dppc.addH.xx
run 200000
Setting up run ...
Memory usage per processor = 105.505 Mbytes
Step CPU Temp Press pull[3] pull[4] pull[5] pull[6] pull[7]
       0 0 323.2253 1.0353196e+15 3.5527137e-13
-3.5527137e-13 5.8980485 5.8996228 -3.5527137e-16
      10 4.697108 1206.4801 -1.2635928e+15 0.032503132
-0.032503132 5.8880485 5.8889738 -0.00012854972
      20 9.5625679 nan nan -1.0421669e-316
1.0421669e-316 5.8780485 nan -0.00072789367
      30 10.417853 nan nan -1.0421669e-316
1.0421669e-316 5.8680485 nan -0.00072789367
      40 11.277054 nan nan -1.0421669e-316
1.0421669e-316 5.8580485 nan -0.00072789367
      50 12.133108 nan nan -1.0421669e-316
1.0421669e-316 5.8480485 nan -0.00072789367
      60 12.990315 nan nan -1.0421669e-316
1.0421669e-316 5.8380485 nan -0.00072789367
      70 13.847533 nan nan -1.0421669e-316
1.0421669e-316 5.8280485 nan -0.00072789367
      80 14.706971 nan nan -1.0421669e-316
1.0421669e-316 5.8180485 nan -0.00072789367
      90 15.563443 nan nan -1.0421669e-316
1.0421669e-316 5.8080485 nan -0.00072789367
     100 16.418641 nan nan -1.0421669e-316
1.0421669e-316 5.7980485 nan -0.00072789367

Di Cheng

University of Science and Technology of China
Hefei, Anhui Province 230026
P. R. China
E-mail: [email protected]...
Tel.: +86-15321055911

3

Hi, Axel

Now I use constant velocity mode, set velocity to be 1.0. Which end of

a velocity of 1.0 is _far_ too large. the velocity unit is one
length unit per time unit, i.e. one angstrom per femtosecond.
you want something like that is at least 4 orders of magnitude
smaller.

the spring is moving, sio2 end or ref.point end?

this is not how it works. at the time when fix smd is
defined, it records the distance between the reference
point and the center of mass of your fix group. this
is the initial value of R. then in each step this is
reduced at the rate that you have given.

[...]

# hold lower lipid bilayer
fix hold_lipid lower_lipid spring/self 10
# pull
fix pull sio2 smd cvel 50 -0.001 tether NULL NULL 37 0
# fix 2_sio2 sio2 nve
fix sio2 sio2 setforce 0 0 NULL

ok. so you want to pull it in a straight line.
is there a specific reason for that?

velocity sio2 set 0 0 -0.001 units box

this is not really needed.

thermo 10
thermo_style custom step cpu temp press f_pull[3] f_pull[4]
f_pull[5] f_pull[6] f_pull[7]

# 1fs per step
timestep 1

dump 1 all atom 10 dump.dppc.addH.xx
run 200000
Setting up run ...
Memory usage per processor = 105.505 Mbytes
Step CPU Temp Press pull[3] pull[4] pull[5] pull[6] pull[7]
       0 0 323.2253 1.0353196e+15 3.5527137e-13

-3.5527137e-13 5.8980485 5.8996228 -3.5527137e-16
      10 4.697108 1206.4801 -1.2635928e+15 0.032503132

have you equilibrated your system before you start the pulling?
you have a gigantic temperature after 10 steps already.
does the temperature explode in the same way when you use
0.0 as pulling velocity?

-0.032503132 5.8880485 5.8889738 -0.00012854972
      20 9.5625679 nan nan -1.0421669e-316
1.0421669e-316 5.8780485 nan -0.00072789367
      30 10.417853 nan nan -1.0421669e-316

at this point all hope is lost...

4

Dear Axel

Thanks for you reply.

the data file is generated by restart2data from a restart, which is
after minimized and 150ps run.
Some other run based on it is OK.

Hi, Axel

Now I use constant velocity mode, set velocity to be 1.0. Which end of

a velocity of 1.0 is _far_ too large. the velocity unit is one
length unit per time unit, i.e. one angstrom per femtosecond.
you want something like that is at least 4 orders of magnitude
smaller.

I mean for example, the velocity is 1.0, which end is moving with 1.0 velocity.
Actually, I am using a velocity of 0.01 to test the in.script now.

the spring is moving, sio2 end or ref.point end?

this is not how it works. at the time when fix smd is
defined, it records the distance between the reference
point and the center of mass of your fix group. this
is the initial value of R. then in each step this is
reduced at the rate that you have given.

[...]

# hold lower lipid bilayer
fix hold_lipid lower_lipid spring/self 10
# pull
fix pull sio2 smd cvel 50 -0.001 tether NULL NULL 37 0
# fix 2_sio2 sio2 nve
fix sio2 sio2 setforce 0 0 NULL

ok. so you want to pull it in a straight line.
is there a specific reason for that?

I want to simulate a AFM pulling: attach a sio2 sphere to AFM tip, and
move it vertically to the membrane.

velocity sio2 set 0 0 -0.001 units box

this is not really needed.

I want to make the initial speed of sio2 equal to speed of spring...

thermo 10
thermo_style custom step cpu temp press f_pull[3] f_pull[4]
f_pull[5] f_pull[6] f_pull[7]

# 1fs per step
timestep 1

dump 1 all atom 10 dump.dppc.addH.xx
run 200000
Setting up run ...
Memory usage per processor = 105.505 Mbytes
Step CPU Temp Press pull[3] pull[4] pull[5] pull[6] pull[7]
0 0 323.2253 1.0353196e+15 3.5527137e-13

-3.5527137e-13 5.8980485 5.8996228 -3.5527137e-16
10 4.697108 1206.4801 -1.2635928e+15 0.032503132

have you equilibrated your system before you start the pulling?
you have a gigantic temperature after 10 steps already.
does the temperature explode in the same way when you use
0.0 as pulling velocity?

The temperature is not being calculated correctly here because not all
groups are fixed nvt. It must be compute in another way. And I do not
really care about it now. What I really want to deal with is to move
the sio2.

-0.032503132 5.8880485 5.8889738 -0.00012854972
20 9.5625679 nan nan -1.0421669e-316
1.0421669e-316 5.8780485 nan -0.00072789367
30 10.417853 nan nan -1.0421669e-316

at this point all hope is lost...

YES, HOPE is lost for this run. So I modified the in.script and try it again

5

Dear Axel

Thanks for you reply.

the data file is generated by restart2data from a restart, which is
after minimized and 150ps run.
Some other run based on it is OK.

ok.

Actually, I am using a velocity of 0.01 to test the in.script now.

that is still very fast.

fix sio2 sio2 setforce 0 0 NULL

ok. so you want to pull it in a straight line.
is there a specific reason for that?

I want to simulate a AFM pulling: attach a sio2 sphere to AFM tip, and
move it vertically to the membrane.

yes. but how strong is the sphere attached to the tip?
in part this is modelled with the spring constant in fix spring.
but if the sphere is actually moving around, you could
emulate it by setting R0 to a finite value. but then again,
if you don't consider the tip in your model and the shockwave
it generates, then you are not that realistic anyways.
the way how i would do this is by building a proper tip,
set all its forces to zero, its velocity to a constant. integrate it
with nve and then attach the sphere to its tip atom(s) with
fix spring couple.

i am currently having some fun doing something very similar
http://sites.google.com/site/akohlmey/software/vrpn-icms
only i'm trying to get it to work interactively... :wink:

velocity sio2 set 0 0 -0.001 units box

this is not really needed.

I want to make the initial speed of sio2 equal to speed of spring...

are you sure it is in the right direction?
if not, the sphere will snap back at much
higher speed.

[...]

The temperature is not being calculated correctly here because not all
groups are fixed nvt. It must be compute in another way. And I do not
really care about it now. What I really want to deal with is to move
the sio2.

nevertheless, i would have a look at it. setting up the temperature
being computed correctly. will help you to see what is going down
and whether you are pulling to hard.

cheers,
    axel.

6

Hi, Axel

I do not want to model a AFM tip. I think a spring with a force
constant of the AFM cantilever is enough for me to simulate system.

And I move the equilibrium position down, the spring will drag the
sio2 down. Now I work around the fix smd and use fix spring directly.

variable y loop 200
label up

variable pos equal 43-20+x\*0\.1 \#speed just for fun fix pull sio2 spring tether 5000 NULL NULL {pos} 0 #5000
maybe too big...
thermo 10
thermo_style custom step cpu temp press c_sio2[3] f_sio2[3] f_pull[3]
# output forces acted upon sio2 sphere
run 1000
next y
unfix pull
jump in.addH up

The result is : sio2 is moving! Some parameters are to be reconciled.

Di Cheng

University of Science and Technology of China
Hefei, Anhui Province 230026
P. R. China
E-mail: [email protected]...
Tel.: +86-15321055911

7

Hi, Axel

I do not want to model a AFM tip. I think a spring with a force
constant of the AFM cantilever is enough for me to simulate system.

well, it depends on whether you model this in solvent or in vacuum.
in the first case, i believe you should not leave out the tip, as it
produces a shockwave that you would be missing otherwise. in
the second case, it would be ok.

And I move the equilibrium position down, the spring will drag the
sio2 down. Now I work around the fix smd and use fix spring directly.

if it works for you, then fine, but this is _exactly_ what fix smd does,
only that it acts as if you'd be using a step size of 1 instead of 200.

cheers,
    axel.