These are fix styles contributed by users, which can be used if +LAMMPS is built with the appropriate package. +
+| smd + |
Compute styles. See the compute command for one-line diff -Naur lammps-18May08/doc/Section_commands.txt lammps-19May08/doc/Section_commands.txt --- lammps-18May08/doc/Section_commands.txt 2008-05-15 15:21:50.000000000 -0600 +++ lammps-19May08/doc/Section_commands.txt 2008-05-15 16:05:51.000000000 -0600 @@ -437,6 +437,11 @@ "wall/reflect"_fix_wall_reflect.html, "wiggle"_fix_wiggle.html :tb(c=8,ea=c) +These are fix styles contributed by users, which can be used if +"LAMMPS is built with the appropriate package"_Section_start.html#2_3. + +"smd"_fix_smd.html :tb(c=4,ea=c) + :line Compute styles. See the "compute"_compute.html command for one-line diff -Naur lammps-18May08/doc/compute_ackland_atom.html lammps-19May08/doc/compute_ackland_atom.html --- lammps-18May08/doc/compute_ackland_atom.html 2008-01-03 17:56:10.000000000 -0700 +++ lammps-19May08/doc/compute_ackland_atom.html 2008-05-15 16:05:51.000000000 -0600 @@ -56,7 +56,11 @@ this section for an overview of LAMMPS output options.
-Restrictions: none +
Restrictions: +
+This compute is part of the "user-ackland" package. It is only +enabled if LAMMPS was built with that package. See the Making +LAMMPS section for more info.
Related commands:
diff -Naur lammps-18May08/doc/compute_ackland_atom.txt lammps-19May08/doc/compute_ackland_atom.txt --- lammps-18May08/doc/compute_ackland_atom.txt 2008-01-03 17:56:10.000000000 -0700 +++ lammps-19May08/doc/compute_ackland_atom.txt 2008-05-15 16:05:51.000000000 -0600 @@ -53,7 +53,11 @@ "this section"_Section_howto.html#4_15 for an overview of LAMMPS output options. -[Restrictions:] none +[Restrictions:] + +This compute is part of the "user-ackland" package. It is only +enabled if LAMMPS was built with that package. See the "Making +LAMMPS"_Section_start.html#2_3 section for more info. [Related commands:] diff -Naur lammps-18May08/doc/fix_drag.html lammps-19May08/doc/fix_drag.html --- lammps-18May08/doc/fix_drag.html 2008-01-02 12:25:15.000000000 -0700 +++ lammps-19May08/doc/fix_drag.html 2008-05-15 16:05:51.000000000 -0600 @@ -59,7 +59,8 @@Related commands:
-fix spring, fix spring/self, +fix spring/rg, fix smd
Default: none
diff -Naur lammps-18May08/doc/fix_drag.txt lammps-19May08/doc/fix_drag.txt --- lammps-18May08/doc/fix_drag.txt 2008-01-02 12:25:15.000000000 -0700 +++ lammps-19May08/doc/fix_drag.txt 2008-05-15 16:05:51.000000000 -0600 @@ -57,6 +57,7 @@ [Related commands:] -"fix spring"_fix_spring.html +"fix spring"_fix_spring.html, "fix spring/self"_fix_spring_self.html, +"fix spring/rg"_fix_spring_rg.html, "fix smd"_fix_smd.html [Default:] none diff -Naur lammps-18May08/doc/fix_smd.html lammps-19May08/doc/fix_smd.html --- lammps-18May08/doc/fix_smd.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/doc/fix_smd.html 2008-05-15 16:05:51.000000000 -0600 @@ -0,0 +1,162 @@ + +Syntax: +
+fix ID group-ID smd type values keyword values ++
tether values = x y z R0 + x,y,z = point to which spring is tethered + R0 = distance of end of spring from tether point (distance units) + couple values = group-ID2 x y z R0 + group-ID2 = 2nd group to couple to fix group with a spring + x,y,z = direction of spring, automatically computed with 'auto' + R0 = distance of end of spring (distance units) ++ +
Examples: +
+fix pull cterm smd cvel 20.0 -0.00005 tether NULL NULL 100.0 0.0 +fix pull cterm smd cvel 20.0 -0.0001 tether 25.0 25 25.0 0.0 +fix stretch cterm smd cvel 20.0 0.0001 couple nterm auto auto auto 0.0 +fix pull cterm smd cfor 5.0 tether 25.0 25.0 25.0 0.0 ++
Description: +
+This fix implements several options of steered MD (SMD) as reviewed in +(Izrailev), which allows to induce conformational changes +in systems and to compute the potential of mean force (PMF) along the +assumed reaction coordinate (Park) based on Jarzynski's +equality (Jarzynski). This fix borrows a lot from fix +spring and fix setforce. +
+You can apply a moving spring force to a group of atoms (tether +style) or between two groups of atoms (couple style). The spring +can then be used in either constant velocity (cvel) mode or in +constant force (cfor) mode to induce transitions in your systems. +When running in tether style, you may need some way to fix some +other part of the system (e.g. via fix +spring/self) +
+The tether style attaches a spring between a point at a distance of +R0 away from a fixed point x,y,z and the center of mass of the fix +group of atoms. A restoring force of magnitude K (R - R0) Mi / M is +applied to each atom in the group where K is the spring constant, Mi +is the mass of the atom, and M is the total mass of all atoms in the +group. Note that K thus represents the total force on the group of +atoms, not a per-atom force. +
+In cvel mode the distance R is incremented or decremented +monotonously according to the pulling (or pushing) velocity. +In cfor mode a constant force is added and the actual distance +in direction of the spring is recorded. +
+The couple style links two groups of atoms together. The first +group is the fix group; the second is specified by group-ID2. The +groups are coupled together by a spring that is at equilibrium when +the two groups are displaced by a vector in direction x,y,z with +respect to each other and at a distance R0 from that displacement. +Note that x,y,z only provides a direction and will be internally +normalized. But since it represents the absolute displacement of +group-ID2 relative to the fix group, (1,1,0) is a different spring +than (-1,-1,0). For each vector component, the displacement can be +described with the auto parameter. In this case the direction is +recomputed in every step, which can be useful for steering a local +process where the whole object undergoes some other change. When the +relative positions and distance between the two groups are not in +equilibrium, the same spring force described above is applied to atoms +in each of the two groups. +
+For both the tether and couple styles, any of the x,y,z values can +be specified as NULL which means do not include that dimension in the +distance calculation or force application. +
+For constant velocity pulling (cvel mode), the running integral +over the pulling force in direction of the spring is recorded and +can then later be used to compute the potential of mean force (PMF) +by averaging over multiple independent trajectories along the same +pulling path. +
+Restart, fix_modify, output, run start/stop, minimize info: +
+The fix stores the direction of the spring, current pulling target +distance and the running PMF to binary restart files. +See the read_restart command for info on how to +re-specify a fix in an input script that reads a restart file, so that +the operation of the fix continues in an uninterrupted fashion. +
+None of the fix_modify options are relevant to this +fix. +
+This fix computes a vector list of 7 quantities, which can be accessed +by various output commands. The quantities +in the vector are in this order: the x-, y-, and z-component of the +pulling force, the total force in direction of the pull, the +equilibrium distance of the spring, the distance between the two +reference points, and finally the accumulated PMF (the sum of pulling +forces times displacement). +
+The force is the total force on the group of atoms by the spring. In +the case of the couple style, it is the force on the fix group +(group-ID) or the negative of the force on the 2nd group (group-ID2). +The vector values calculated by this fix are "extensive", meaning they +scale with the number of atoms in the simulation. +
+No parameter of this fix can be used with the start/stop keywords of +the run command. This fix is not invoked during energy +minimization. +
+Restrictions: +
+This fix is part of the "user-smd" package. It is only enabled if +LAMMPS was built with that package. See the Making +LAMMPS section for more info. +
+Related commands: +
+fix drag, fix spring, +fix spring/self, +fix spring/rg +
+Default: none +
+(Izrailev) Izrailev, Stepaniants, Isralewitz, Kosztin, Lu, Molnar, +Wriggers, Schulten. Computational Molecular Dynamics: Challenges, +Methods, Ideas, volume 4 of Lecture Notes in Computational Science and +Engineering, pp. 39-65. Springer-Verlag, Berlin, 1998. +
+(Park) +Park, Schulten, J. Chem. Phys. 120 (13), 5946 (2004) +
+(Jarzynski) +Jarzynski, Phys. Rev. Lett. 78, 2690 (1997) +
+ diff -Naur lammps-18May08/doc/fix_smd.txt lammps-19May08/doc/fix_smd.txt --- lammps-18May08/doc/fix_smd.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/doc/fix_smd.txt 2008-05-15 16:05:51.000000000 -0600 @@ -0,0 +1,151 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix smd command :h3 + +[Syntax:] + +fix ID group-ID smd type values keyword values :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +smd = style name of this fix command :l +mode = {cvel} or {cfor} to select constant velocity or constant force SMD :l + {cvel} values = K vel + K = spring constant (force/distance units) + vel = velocity of pulling (distance/time units) + {cfor} values = force + force = pulling force (force units) +keyword = {tether} or {couple} :l + {tether} values = x y z R0 + x,y,z = point to which spring is tethered + R0 = distance of end of spring from tether point (distance units) + {couple} values = group-ID2 x y z R0 + group-ID2 = 2nd group to couple to fix group with a spring + x,y,z = direction of spring, automatically computed with 'auto' + R0 = distance of end of spring (distance units) :pre +:ule + +[Examples:] + +fix pull cterm smd cvel 20.0 -0.00005 tether NULL NULL 100.0 0.0 +fix pull cterm smd cvel 20.0 -0.0001 tether 25.0 25 25.0 0.0 +fix stretch cterm smd cvel 20.0 0.0001 couple nterm auto auto auto 0.0 +fix pull cterm smd cfor 5.0 tether 25.0 25.0 25.0 0.0 :pre + +[Description:] + +This fix implements several options of steered MD (SMD) as reviewed in +"(Izrailev)"_#Izrailev, which allows to induce conformational changes +in systems and to compute the potential of mean force (PMF) along the +assumed reaction coordinate "(Park)"_#Park based on Jarzynski's +equality "(Jarzynski)"_#Jarzynski. This fix borrows a lot from "fix +spring"_fix_spring.html and "fix setforce"_fix_setforce.html. + +You can apply a moving spring force to a group of atoms ({tether} +style) or between two groups of atoms ({couple} style). The spring +can then be used in either constant velocity ({cvel}) mode or in +constant force ({cfor}) mode to induce transitions in your systems. +When running in {tether} style, you may need some way to fix some +other part of the system (e.g. via "fix +spring/self"_fix_spring_self.html) + +The {tether} style attaches a spring between a point at a distance of +R0 away from a fixed point {x,y,z} and the center of mass of the fix +group of atoms. A restoring force of magnitude K (R - R0) Mi / M is +applied to each atom in the group where {K} is the spring constant, Mi +is the mass of the atom, and M is the total mass of all atoms in the +group. Note that {K} thus represents the total force on the group of +atoms, not a per-atom force. + +In {cvel} mode the distance R is incremented or decremented +monotonously according to the pulling (or pushing) velocity. +In {cfor} mode a constant force is added and the actual distance +in direction of the spring is recorded. + +The {couple} style links two groups of atoms together. The first +group is the fix group; the second is specified by group-ID2. The +groups are coupled together by a spring that is at equilibrium when +the two groups are displaced by a vector in direction {x,y,z} with +respect to each other and at a distance R0 from that displacement. +Note that {x,y,z} only provides a direction and will be internally +normalized. But since it represents the {absolute} displacement of +group-ID2 relative to the fix group, (1,1,0) is a different spring +than (-1,-1,0). For each vector component, the displacement can be +described with the {auto} parameter. In this case the direction is +recomputed in every step, which can be useful for steering a local +process where the whole object undergoes some other change. When the +relative positions and distance between the two groups are not in +equilibrium, the same spring force described above is applied to atoms +in each of the two groups. + +For both the {tether} and {couple} styles, any of the x,y,z values can +be specified as NULL which means do not include that dimension in the +distance calculation or force application. + +For constant velocity pulling ({cvel} mode), the running integral +over the pulling force in direction of the spring is recorded and +can then later be used to compute the potential of mean force (PMF) +by averaging over multiple independent trajectories along the same +pulling path. + +[Restart, fix_modify, output, run start/stop, minimize info:] + +The fix stores the direction of the spring, current pulling target +distance and the running PMF to "binary restart files"_restart.html. +See the "read_restart"_read_restart.html command for info on how to +re-specify a fix in an input script that reads a restart file, so that +the operation of the fix continues in an uninterrupted fashion. + +None of the "fix_modify"_fix_modify.html options are relevant to this +fix. + +This fix computes a vector list of 7 quantities, which can be accessed +by various "output commands"_Section_howto.html#4_15. The quantities +in the vector are in this order: the x-, y-, and z-component of the +pulling force, the total force in direction of the pull, the +equilibrium distance of the spring, the distance between the two +reference points, and finally the accumulated PMF (the sum of pulling +forces times displacement). + +The force is the total force on the group of atoms by the spring. In +the case of the {couple} style, it is the force on the fix group +(group-ID) or the negative of the force on the 2nd group (group-ID2). +The vector values calculated by this fix are "extensive", meaning they +scale with the number of atoms in the simulation. + +No parameter of this fix can be used with the {start/stop} keywords of +the "run"_run.html command. This fix is not invoked during "energy +minimization"_minimize.html. + +[Restrictions:] + +This fix is part of the "user-smd" package. It is only enabled if +LAMMPS was built with that package. See the "Making +LAMMPS"_Section_start.html#2_3 section for more info. + +[Related commands:] + +"fix drag"_fix_drag.html, "fix spring"_fix_spring.html, +"fix spring/self"_fix_spring_self.html, +"fix spring/rg"_fix_spring_rg.html + +[Default:] none + +:line + +:link(Israilev) +[(Izrailev)] Izrailev, Stepaniants, Isralewitz, Kosztin, Lu, Molnar, +Wriggers, Schulten. Computational Molecular Dynamics: Challenges, +Methods, Ideas, volume 4 of Lecture Notes in Computational Science and +Engineering, pp. 39-65. Springer-Verlag, Berlin, 1998. + +[(Park)] +Park, Schulten, J. Chem. Phys. 120 (13), 5946 (2004) + +[(Jarzynski)] +Jarzynski, Phys. Rev. Lett. 78, 2690 (1997) diff -Naur lammps-18May08/doc/fix_spring.html lammps-19May08/doc/fix_spring.html --- lammps-18May08/doc/fix_spring.html 2008-01-02 12:25:15.000000000 -0700 +++ lammps-19May08/doc/fix_spring.html 2008-05-15 16:05:51.000000000 -0600 @@ -113,8 +113,8 @@Related commands:
-fix drag, fix spring/self, fix -spring/rg +
fix drag, fix spring/self, +fix spring/rg, fix smd
Default: none
diff -Naur lammps-18May08/doc/fix_spring.txt lammps-19May08/doc/fix_spring.txt --- lammps-18May08/doc/fix_spring.txt 2008-01-02 12:25:15.000000000 -0700 +++ lammps-19May08/doc/fix_spring.txt 2008-05-15 16:05:51.000000000 -0600 @@ -106,7 +106,8 @@ [Related commands:] -"fix drag"_fix_drag.html, "fix spring/self"_fix_spring_self.html, "fix -spring/rg"_fix_spring_rg.html +"fix drag"_fix_drag.html, "fix spring/self"_fix_spring_self.html, +"fix spring/rg"_fix_spring_rg.html, "fix smd"_fix_smd.html + [Default:] none diff -Naur lammps-18May08/doc/fix_spring_rg.html lammps-19May08/doc/fix_spring_rg.html --- lammps-18May08/doc/fix_spring_rg.html 2007-10-10 16:28:11.000000000 -0600 +++ lammps-19May08/doc/fix_spring_rg.html 2008-05-15 16:05:51.000000000 -0600 @@ -65,6 +65,7 @@Related commands:
fix spring, fix spring/self +fix drag, fix smd
Default: none
diff -Naur lammps-18May08/doc/fix_spring_rg.txt lammps-19May08/doc/fix_spring_rg.txt --- lammps-18May08/doc/fix_spring_rg.txt 2007-10-10 16:28:11.000000000 -0600 +++ lammps-19May08/doc/fix_spring_rg.txt 2008-05-15 16:05:51.000000000 -0600 @@ -61,5 +61,6 @@ [Related commands:] "fix spring"_fix_spring.html, "fix spring/self"_fix_spring_self.html +"fix drag"_fix_drag.html, "fix smd"_fix_smd.html [Default:] none diff -Naur lammps-18May08/doc/fix_spring_self.html lammps-19May08/doc/fix_spring_self.html --- lammps-18May08/doc/fix_spring_self.html 2007-10-10 16:28:11.000000000 -0600 +++ lammps-19May08/doc/fix_spring_self.html 2008-05-15 16:05:51.000000000 -0600 @@ -52,7 +52,8 @@Related commands:
-fix drag, fix spring, +fix smd, fix spring/rg
Default: none
diff -Naur lammps-18May08/doc/fix_spring_self.txt lammps-19May08/doc/fix_spring_self.txt --- lammps-18May08/doc/fix_spring_self.txt 2007-10-10 16:28:11.000000000 -0600 +++ lammps-19May08/doc/fix_spring_self.txt 2008-05-15 16:05:51.000000000 -0600 @@ -49,6 +49,7 @@ [Related commands:] -"fix drag"_fix_drag.html, "fix spring"_fix_spring.html +"fix drag"_fix_drag.html, "fix spring"_fix_spring.html, +"fix smd"_fix_smd.html, "fix spring/rg"_fix_spring_rg.html [Default:] none diff -Naur lammps-18May08/src/Makefile lammps-19May08/src/Makefile --- lammps-18May08/src/Makefile 2008-02-22 07:28:58.000000000 -0700 +++ lammps-19May08/src/Makefile 2008-05-15 16:06:17.000000000 -0600 @@ -16,7 +16,7 @@ PACKAGE = asphere class2 colloid dipole dpd granular \ kspace manybody meam molecule opt poems xtc -PACKUSER = user-ackland user-cg-cmm user-ewaldn +PACKUSER = user-ackland user-cg-cmm user-ewaldn user-smd PACKALL = $(PACKAGE) $(PACKUSER) diff -Naur lammps-18May08/src/USER-SMD/Install.csh lammps-19May08/src/USER-SMD/Install.csh --- lammps-18May08/src/USER-SMD/Install.csh 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/src/USER-SMD/Install.csh 2008-05-15 16:06:17.000000000 -0600 @@ -0,0 +1,20 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp -p style_user_smd.h .. + + cp -p fix_smd.cpp .. + + cp -p fix_smd.h .. + +else if ($1 == 0) then + + rm ../style_user_smd.h + touch ../style_user_smd.h + + rm ../fix_smd.cpp + + rm ../fix_smd.h + +endif diff -Naur lammps-18May08/src/USER-SMD/README lammps-19May08/src/USER-SMD/README --- lammps-18May08/src/USER-SMD/README 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/src/USER-SMD/README 2008-05-15 16:08:14.000000000 -0600 @@ -0,0 +1,20 @@ +The files in this directory are a user-contributed package for LAMMPS. + +The person who created these files is Axel Kohlmeyer +(akohlmey@cmm.chem.upenn.edu). Contact him directly if you have +questions or for example scripts that use it. + +This package implements a "fix smd" command which can be used in a +LAMMPS input script. SMD stands for steered molecular dynamics. + +Here is what Axel said about the new fix: + +Please find attached a fix that I was asked to write for our +"coarse-grainers". They need free energy profiles for their +parametrization and this seems to work best. The code basically is a +version of "fix spring" on steroids (constant velocity pulling with +velocity set to zero should recover "fix spring"), but I found it +cleaner to have it in a seprarate class since otherwise the result may +become too convoluted for people that only want to use "fix spring". + + diff -Naur lammps-18May08/src/USER-SMD/fix_smd.cpp lammps-19May08/src/USER-SMD/fix_smd.cpp --- lammps-18May08/src/USER-SMD/fix_smd.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/src/USER-SMD/fix_smd.cpp 2008-05-15 16:06:17.000000000 -0600 @@ -0,0 +1,406 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (UPenn) + based on fix spring by: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_smd.h" +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "error.h" +#include "group.h" + +using namespace LAMMPS_NS; + +enum { SMD_NONE=0, + SMD_TETHER=1<<0, SMD_COUPLE=1<<1, + SMD_CVEL=1<<2, SMD_CFOR=1<<3, + SMD_AUTOX=1<<4, SMD_AUTOY=1<<5, SMD_AUTOZ=1<<6}; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + styleflag = SMD_NONE; + k_smd = f_smd = v_smd = -1.0; + xflag = yflag = zflag = 1; + xc = yc = zc = 0.0; + xn = yn = zn = 1.0; + pmf = r_old = r_now = r0 = 0.0; + + restart_global = 1; + vector_flag = 1; + size_vector = 7; + scalar_vector_freq = 1; + extvector = 1; + + int argoffs=3; + if (strcmp(arg[argoffs],"cvel") == 0) { + if (narg < argoffs+3) error->all("Illegal fix smd command"); + styleflag |= SMD_CVEL; + k_smd = atof(arg[argoffs+1]); + v_smd = atof(arg[argoffs+2]); // to be multiplied by update->dt when used. + argoffs += 3; + } else if (strcmp(arg[argoffs],"cfor") == 0) { + if (narg < argoffs+2) error->all("Illegal fix smd command"); + styleflag |= SMD_CFOR; + f_smd = atof(arg[argoffs+1]); + argoffs += 2; + } else error->all("Illegal fix smd command"); + + if (strcmp(arg[argoffs],"tether") == 0) { + if (narg < argoffs+5) error->all("Illegal fix smd command"); + styleflag |= SMD_TETHER; + if (strcmp(arg[argoffs+1],"NULL") == 0) xflag = 0; + else xc = atof(arg[argoffs+1]); + if (strcmp(arg[argoffs+2],"NULL") == 0) yflag = 0; + else yc = atof(arg[argoffs+2]); + if (strcmp(arg[argoffs+3],"NULL") == 0) zflag = 0; + else zc = atof(arg[argoffs+3]); + r0 = atof(arg[argoffs+4]); + if (r0 < 0) error->all("R0 < 0 for fix smd command"); + argoffs += 5; + } else if (strcmp(arg[argoffs],"couple") == 0) { + if (narg < argoffs+6) error->all("Illegal fix smd command"); + styleflag |= SMD_COUPLE; + igroup2 = group->find(arg[argoffs+1]); + if (igroup2 == -1) + error->all("Could not find fix smd couple group ID"); + if (igroup2 == igroup) + error->all("Two groups cannot be the same in fix smd couple"); + group2bit = group->bitmask[igroup2]; + + if (strcmp(arg[argoffs+2],"NULL") == 0) xflag = 0; + else if (strcmp(arg[argoffs+2],"auto") == 0) styleflag |= SMD_AUTOX; + else xc = atof(arg[argoffs+2]); + if (strcmp(arg[argoffs+3],"NULL") == 0) yflag = 0; + else if (strcmp(arg[argoffs+3],"auto") == 0) styleflag |= SMD_AUTOY; + else yc = atof(arg[argoffs+3]); + if (strcmp(arg[argoffs+4],"NULL") == 0) zflag = 0; + else if (strcmp(arg[argoffs+4],"auto") == 0) styleflag |= SMD_AUTOZ; + else zc = atof(arg[argoffs+4]); + + r0 = atof(arg[argoffs+5]); + if (r0 < 0) error->all("R0 < 0 for fix smd command"); + argoffs +=6; + } else error->all("Illegal fix smd command"); + + force_flag = 0; + ftotal[0] = ftotal[1] = ftotal[2] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int FixSMD::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::init() +{ + double xcm[3], xcm2[3]; + masstotal = group->mass(igroup); + group->xcm(igroup,masstotal,xcm); + + double dx,dy,dz; + if (styleflag & SMD_TETHER) { + dx = xc - xcm[0]; + dy = yc - xcm[1]; + dz = zc - xcm[2]; + } else { /* SMD_COUPLE */ + masstotal2 = group->mass(igroup2); + group->xcm(igroup2,masstotal2,xcm2); + if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0]; + else dx = xc; + if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1]; + else dy = yc; + if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2]; + else dz = zc; + } + + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r_old = sqrt(dx*dx + dy*dy + dz*dz); + if (r_old > SMALL) { + xn = dx/r_old; + yn = dy/r_old; + zn = dz/r_old; + } + + if (strcmp(update->integrate_style,"respa") == 0) + nlevels_respa = ((Respa *) update->integrate)->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::setup(int vflag) +{ + if (strcmp(update->integrate_style,"verlet") == 0) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::post_force(int vflag) +{ + if (styleflag & SMD_TETHER) smd_tether(); + else smd_couple(); + + if (styleflag & SMD_CVEL) r_old += v_smd * update->dt; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::smd_tether() +{ + double xcm[3]; + group->xcm(igroup,masstotal,xcm); + + // fx,fy,fz = components of k * (r-r0) + + double dx,dy,dz,fx,fy,fz,r,dr; + + dx = xcm[0] - xc; + dy = xcm[1] - yc; + dz = xcm[2] - zc; + r_now = sqrt(dx*dx + dy*dy + dz*dz); + + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r = sqrt(dx*dx + dy*dy + dz*dz); + if (styleflag & SMD_CVEL) { + if(r > SMALL) { + double fsign; + fsign = (v_smd<0.0) ? -1.0 : 1.0; + dr = r - r0 - r_old; + fx = k_smd*dx*dr/r; + fy = k_smd*dy*dr/r; + fz = k_smd*dz*dr/r; + pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt; + } + } else { + r_old = r; + fx = f_smd*dx/r; + fy = f_smd*dy/r; + fz = f_smd*dz/r; + } + + // apply restoring force to atoms in group + // f = -k*(r-r0)*mass/masstotal + + double **f = atom->f; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + int nlocal = atom->nlocal; + + ftotal[0] = ftotal[1] = ftotal[2] = 0.0; + force_flag = 0; + + double massfrac; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + massfrac = mass[type[i]]/masstotal; + fx *= massfrac; + fy *= massfrac; + fz *= massfrac; + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + ftotal[0] -= fx; + ftotal[1] -= fy; + ftotal[2] -= fz; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::smd_couple() +{ + double xcm[3],xcm2[3]; + group->xcm(igroup,masstotal,xcm); + group->xcm(igroup2,masstotal2,xcm2); + + // renormalize direction of spring + double dx,dy,dz,r,dr; + if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0]; + else dx = xn*r_old; + if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1]; + else dy = yn*r_old; + if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2]; + else dz = zn*r_old; + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r = sqrt(dx*dx + dy*dy + dz*dz); + if (r > SMALL) { + xn = dx/r; yn = dy/r; zn = dz/r; + } + + double fx,fy,fz; + if (styleflag & SMD_CVEL) { + dx = xcm2[0] - xcm[0]; + dy = xcm2[1] - xcm[1]; + dz = xcm2[2] - xcm[2]; + r_now = sqrt(dx*dx + dy*dy + dz*dz); + + dx -= xn*r_old; + dy -= yn*r_old; + dz -= zn*r_old; + + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r = sqrt(dx*dx + dy*dy + dz*dz); + dr = r - r0; + + if (r > SMALL) { + double fsign; + fsign = (v_smd<0.0) ? -1.0 : 1.0; + + fx = k_smd*dx*dr/r; + fy = k_smd*dy*dr/r; + fz = k_smd*dz*dr/r; + pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt; + } + + } else { + dx = xcm2[0] - xcm[0]; + dy = xcm2[1] - xcm[1]; + dz = xcm2[2] - xcm[2]; + r_now = sqrt(dx*dx + dy*dy + dz*dz); + r_old = r; + + fx = f_smd*xn; + fy = f_smd*yn; + fz = f_smd*zn; + } + + // apply restoring force to atoms in group + // f = -k*(r-r0)*mass/masstotal + + double **f = atom->f; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + int nlocal = atom->nlocal; + + ftotal[0] = ftotal[1] = ftotal[2] = 0.0; + force_flag = 0; + + double massfrac; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + massfrac = mass[type[i]]/masstotal; + f[i][0] += fx*massfrac; + f[i][1] += fy*massfrac; + f[i][2] += fz*massfrac; + ftotal[0] += fx*massfrac; + ftotal[1] += fy*massfrac; + ftotal[2] += fz*massfrac; + } + if (mask[i] & group2bit) { + massfrac = mass[type[i]]/masstotal2; + f[i][0] -= fx*massfrac; + f[i][1] -= fy*massfrac; + f[i][2] -= fz*massfrac; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::write_restart(FILE *fp) +{ +#define RESTART_ITEMS 5 + double buf[RESTART_ITEMS], fsign; + + if (comm->me == 0) { + // make sure we project the force into the direction of the pulling. + fsign = (v_smd<0.0) ? -1.0 : 1.0; + buf[0] = r_old; + buf[1] = xn*fsign; + buf[2] = yn*fsign; + buf[3] = zn*fsign; + buf[4] = pmf; + int size = RESTART_ITEMS*sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(&buf[0],sizeof(double),RESTART_ITEMS,fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::restart(char *buf) +{ + double *list = (double *)buf; + r_old = list[0]; + xn=list[1]; + yn=list[2]; + zn=list[3]; + pmf=list[4]; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- + return components of total smd force on fix group +------------------------------------------------------------------------- */ + +double FixSMD::compute_vector(int n) +{ + // only sum across procs one time + + if (force_flag == 0) { + MPI_Allreduce(ftotal,ftotal_all,3,MPI_DOUBLE,MPI_SUM,world); + force_flag = 1; + if (styleflag & SMD_CVEL) { + ftotal_all[3]=ftotal_all[0]*xn+ftotal_all[1]*yn+ftotal_all[2]*zn; + ftotal_all[4]=r_old; + } else { + ftotal_all[3]=f_smd; + ftotal_all[4]=r_old; + } + ftotal_all[5]=r_now; + ftotal_all[6]=pmf; + } + return ftotal_all[n]; +} diff -Naur lammps-18May08/src/USER-SMD/fix_smd.h lammps-19May08/src/USER-SMD/fix_smd.h --- lammps-18May08/src/USER-SMD/fix_smd.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/src/USER-SMD/fix_smd.h 2008-05-15 16:06:17.000000000 -0600 @@ -0,0 +1,53 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef FIX_SMD_H +#define FIX_SMD_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixSMD : public Fix { + public: + FixSMD(class LAMMPS *, int, char **); + int setmask(); + void init(); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double compute_vector(int); + + void write_restart(FILE *); + void restart(char *); + + private: + double xc,yc,zc,xn,yn,zn,r0; + double k_smd,f_smd,v_smd; + int xflag,yflag,zflag; + int styleflag; + double r_old,r_now,pmf; + + int igroup2,group2bit; + double masstotal,masstotal2; + int nlevels_respa; + double ftotal[3],ftotal_all[7]; + int force_flag; + + void smd_tether(); + void smd_couple(); +}; + +} + +#endif diff -Naur lammps-18May08/src/USER-SMD/style_user_smd.h lammps-19May08/src/USER-SMD/style_user_smd.h --- lammps-18May08/src/USER-SMD/style_user_smd.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-19May08/src/USER-SMD/style_user_smd.h 2008-05-15 16:06:17.000000000 -0600 @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FixInclude +#include "fix_smd.h" +#endif + +#ifdef FixClass +FixStyle(smd,FixSMD) +#endif diff -Naur lammps-18May08/src/style_user_packages.h lammps-19May08/src/style_user_packages.h --- lammps-18May08/src/style_user_packages.h 2008-02-08 14:57:06.000000000 -0700 +++ lammps-19May08/src/style_user_packages.h 2008-05-15 16:06:17.000000000 -0600 @@ -17,3 +17,4 @@ #include "style_user_ackland.h" #include "style_user_cg_cmm.h" #include "style_user_ewaldn.h" +#include "style_user_smd.h"