diff -Naur lammps-11Nov09/doc/Section_modify.html lammps-12Nov09/doc/Section_modify.html --- lammps-11Nov09/doc/Section_modify.html 2009-06-30 13:06:26.000000000 -0600 +++ lammps-12Nov09/doc/Section_modify.html 2009-11-09 11:37:37.000000000 -0700 @@ -164,11 +164,21 @@ grow re-allocate atom arrays to longer lengths copy copy info for one atom to another atom's array locations pack_comm store an atom's info in a buffer communicated every timestep +pack_comm_vel add velocity info to buffer +pack_comm_one store extra info unique to this atom style unpack_comm retrieve an atom's info from the buffer +unpack_comm_vel also retrieve velocity info +unpack_comm_one retreive extra info unique to this atom style pack_reverse store an atom's info in a buffer communicating partial forces +pack_reverse_one store extra info unique to this atom style unpack_reverse retrieve an atom's info from the buffer +unpack_reverse_one retreive extra info unique to this atom style pack_border store an atom's info in a buffer communicated on neighbor re-builds +pack_border_vel add velocity info to buffer +pack_border_one store extra info unique to this atom style unpack_border retrieve an atom's info from the buffer +unpack_border_vel also retrieve velocity info +unpack_border_one retreive extra info unique to this atom style pack_exchange store all an atom's info to migrate to another processor unpack_exchange retrieve an atom's info from the buffer size_restart number of restart quantities associated with proc's atoms @@ -180,9 +190,10 @@

The constructor of the derived class sets values for several variables -that you must set when defining a new atom style. The atom arrays -themselves are defined in atom.cpp. Search for the word "customize" -and you will find locations you will need to modify. +that you must set when defining a new atom style, which are documented +in atom_vec.h. New atom arrays are defined in atom.cpp. Search for +the word "customize" and you will find locations you will need to +modify.


diff -Naur lammps-11Nov09/doc/Section_modify.txt lammps-12Nov09/doc/Section_modify.txt --- lammps-11Nov09/doc/Section_modify.txt 2009-06-30 13:06:26.000000000 -0600 +++ lammps-12Nov09/doc/Section_modify.txt 2009-11-09 11:37:37.000000000 -0700 @@ -161,11 +161,21 @@ grow: re-allocate atom arrays to longer lengths copy: copy info for one atom to another atom's array locations pack_comm: store an atom's info in a buffer communicated every timestep +pack_comm_vel: add velocity info to buffer +pack_comm_one: store extra info unique to this atom style unpack_comm: retrieve an atom's info from the buffer +unpack_comm_vel: also retrieve velocity info +unpack_comm_one: retreive extra info unique to this atom style pack_reverse: store an atom's info in a buffer communicating partial forces +pack_reverse_one: store extra info unique to this atom style unpack_reverse: retrieve an atom's info from the buffer +unpack_reverse_one: retreive extra info unique to this atom style pack_border: store an atom's info in a buffer communicated on neighbor re-builds +pack_border_vel: add velocity info to buffer +pack_border_one: store extra info unique to this atom style unpack_border: retrieve an atom's info from the buffer +unpack_border_vel: also retrieve velocity info +unpack_border_one: retreive extra info unique to this atom style pack_exchange: store all an atom's info to migrate to another processor unpack_exchange: retrieve an atom's info from the buffer size_restart: number of restart quantities associated with proc's atoms @@ -176,9 +186,10 @@ memory_usage: tally memory allocated by atom arrays :tb(s=:) The constructor of the derived class sets values for several variables -that you must set when defining a new atom style. The atom arrays -themselves are defined in atom.cpp. Search for the word "customize" -and you will find locations you will need to modify. +that you must set when defining a new atom style, which are documented +in atom_vec.h. New atom arrays are defined in atom.cpp. Search for +the word "customize" and you will find locations you will need to +modify. :line diff -Naur lammps-11Nov09/doc/Section_start.html lammps-12Nov09/doc/Section_start.html --- lammps-11Nov09/doc/Section_start.html 2009-11-06 10:12:08.000000000 -0700 +++ lammps-12Nov09/doc/Section_start.html 2009-11-09 11:32:16.000000000 -0700 @@ -360,7 +360,6 @@ class2 class 2 force fields colloid colloidal particle force fields dipole point dipole particles and force fields -dpd dissipative particle dynamics (DPD) force field dsmc Direct Simulation Monte Carlo (DMSC) pair style gpu GPU-enabled force field styles granular force fields and boundary conditions for granular systems diff -Naur lammps-11Nov09/doc/Section_start.txt lammps-12Nov09/doc/Section_start.txt --- lammps-11Nov09/doc/Section_start.txt 2009-11-06 10:12:08.000000000 -0700 +++ lammps-12Nov09/doc/Section_start.txt 2009-11-09 11:32:16.000000000 -0700 @@ -354,7 +354,6 @@ class2 : class 2 force fields colloid : colloidal particle force fields dipole : point dipole particles and force fields -dpd : dissipative particle dynamics (DPD) force field dsmc : Direct Simulation Monte Carlo (DMSC) pair style gpu : GPU-enabled force field styles granular : force fields and boundary conditions for granular systems diff -Naur lammps-11Nov09/doc/atom_style.html lammps-12Nov09/doc/atom_style.html --- lammps-11Nov09/doc/atom_style.html 2009-11-06 14:01:28.000000000 -0700 +++ lammps-12Nov09/doc/atom_style.html 2009-11-09 11:30:55.000000000 -0700 @@ -15,7 +15,7 @@

atom_style style args 
 
- @@ -30,7 +31,8 @@

communicate multi
 communicate multi group solvent
-communicate single cutoff 5.0 
+communicate single ghost yes
+communicate single cutoff 5.0 ghost yes 
 

Description:

@@ -52,13 +54,6 @@ neighbor list construction option that may also be beneficial for simulations of this kind.

-

The group option will limit communication to atoms in the specified -group. This can be useful for models where no ghost atoms are needed -for some kinds of particles. All atoms (not just those in the -specified group) will still migrate to new processors as they move. -The group specified with this option must also be specified via the -atom_modify first command. -

The cutoff option allows you to set a ghost cutoff distance, which is the distance from the borders of a processor's sub-domain at which ghost atoms are acquired from other processors. By default the ghost @@ -98,6 +93,22 @@ potential). Setting the ghost cutoff appropriately can insure it will find the needed atoms.

+

The group option will limit communication to atoms in the specified +group. This can be useful for models where no ghost atoms are needed +for some kinds of particles. All atoms (not just those in the +specified group) will still migrate to new processors as they move. +The group specified with this option must also be specified via the +atom_modify first command. +

+

The vel option enables velocity information to be communicated with +ghost particles. Depending on the atom_style, +velocity info includes the translational velocity, angular velocity, +and angular momentum of a particle. If the vel option is set to +yes, then ghost atoms store these quantities; if no then they do +not. The yes setting is needed by some pair styles which require +the velocity state of both the I and J particles to compute a pairwise +I,J interaction. +

Restrictions: none

Related commands: @@ -106,8 +117,8 @@

Default:

-

The default settings are style = single, group = all, cutoff = 0.0. -The cutoff default of 0.0 means that effectively ghost cutoff = +

The default settings are style = single, group = all, cutoff = 0.0, +ghost = no. The cutoff default of 0.0 means that ghost cutoff = neighbor cutoff = pairwise force cutoff + neighbor skin.

diff -Naur lammps-11Nov09/doc/communicate.txt lammps-12Nov09/doc/communicate.txt --- lammps-11Nov09/doc/communicate.txt 2008-08-21 08:03:21.000000000 -0600 +++ lammps-12Nov09/doc/communicate.txt 2009-11-09 11:30:55.000000000 -0700 @@ -14,16 +14,18 @@ style = {single} or {multi} :ulb,l zero or more keyword/value pairs may be appended :l -keyword = {group} or {cutoff} :l +keyword = {cutoff} or {group} or {vel} :l + {cutoff} value = Rcut (distance units) = communicate atoms from this far away {group} value = group-ID = only communicate atoms in the group - {cutoff} value = Rcut (distance units) = communicate atoms from this far away :pre + {vel} value = {yes} or {no} = do or do not communicate velocity info with ghost atoms :pre :ule [Examples:] communicate multi communicate multi group solvent -communicate single cutoff 5.0 :pre +communicate single ghost yes +communicate single cutoff 5.0 ghost yes :pre [Description:] @@ -45,13 +47,6 @@ neighbor list construction option that may also be beneficial for simulations of this kind. -The {group} option will limit communication to atoms in the specified -group. This can be useful for models where no ghost atoms are needed -for some kinds of particles. All atoms (not just those in the -specified group) will still migrate to new processors as they move. -The group specified with this option must also be specified via the -"atom_modify first"_atom_modify.html command. - The {cutoff} option allows you to set a ghost cutoff distance, which is the distance from the borders of a processor's sub-domain at which ghost atoms are acquired from other processors. By default the ghost @@ -91,6 +86,22 @@ potential). Setting the ghost cutoff appropriately can insure it will find the needed atoms. +The {group} option will limit communication to atoms in the specified +group. This can be useful for models where no ghost atoms are needed +for some kinds of particles. All atoms (not just those in the +specified group) will still migrate to new processors as they move. +The group specified with this option must also be specified via the +"atom_modify first"_atom_modify.html command. + +The {vel} option enables velocity information to be communicated with +ghost particles. Depending on the "atom_style"_atom_style.html, +velocity info includes the translational velocity, angular velocity, +and angular momentum of a particle. If the {vel} option is set to +{yes}, then ghost atoms store these quantities; if {no} then they do +not. The {yes} setting is needed by some pair styles which require +the velocity state of both the I and J particles to compute a pairwise +I,J interaction. + [Restrictions:] none [Related commands:] @@ -99,6 +110,6 @@ [Default:] -The default settings are style = single, group = all, cutoff = 0.0. -The cutoff default of 0.0 means that effectively ghost cutoff = +The default settings are style = single, group = all, cutoff = 0.0, +ghost = no. The cutoff default of 0.0 means that ghost cutoff = neighbor cutoff = pairwise force cutoff + neighbor skin. diff -Naur lammps-11Nov09/doc/pair_dpd.html lammps-12Nov09/doc/pair_dpd.html --- lammps-11Nov09/doc/pair_dpd.html 2009-11-02 11:25:19.000000000 -0700 +++ lammps-12Nov09/doc/pair_dpd.html 2009-11-09 11:32:16.000000000 -0700 @@ -95,10 +95,6 @@

Restrictions:

-

This style is part of the "dpd" package. It is only enabled if -LAMMPS was built with that package. See the Making -LAMMPS section for more info. -

The default frequency for rebuilding neighbor lists is every 10 steps (see the neigh_modify command). This may be too infrequent for DPD simulations since particles move rapidly and can diff -Naur lammps-11Nov09/doc/pair_dpd.txt lammps-12Nov09/doc/pair_dpd.txt --- lammps-11Nov09/doc/pair_dpd.txt 2009-11-02 11:25:19.000000000 -0700 +++ lammps-12Nov09/doc/pair_dpd.txt 2009-11-09 11:32:16.000000000 -0700 @@ -92,10 +92,6 @@ [Restrictions:] -This style is part of the "dpd" 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. - The default frequency for rebuilding neighbor lists is every 10 steps (see the "neigh_modify"_neigh_modify.html command). This may be too infrequent for DPD simulations since particles move rapidly and can diff -Naur lammps-11Nov09/doc/read_data.html lammps-12Nov09/doc/read_data.html --- lammps-11Nov09/doc/read_data.html 2009-11-06 14:01:28.000000000 -0700 +++ lammps-12Nov09/doc/read_data.html 2009-11-09 11:30:55.000000000 -0700 @@ -272,7 +272,6 @@ charge atom-ID atom-type q x y z colloid atom-ID atom-type x y z dipole atom-ID atom-type q x y z mux muy muz -dpd atom-ID atom-type x y z ellipsoid atom-ID atom-type x y z quatw quati quatj quatk full atom-ID molecule-ID atom-type q x y z granular atom-ID atom-type diameter density x y z diff -Naur lammps-11Nov09/doc/read_data.txt lammps-12Nov09/doc/read_data.txt --- lammps-11Nov09/doc/read_data.txt 2009-11-06 14:01:28.000000000 -0700 +++ lammps-12Nov09/doc/read_data.txt 2009-11-09 11:30:55.000000000 -0700 @@ -252,7 +252,6 @@ charge: atom-ID atom-type q x y z colloid: atom-ID atom-type x y z dipole: atom-ID atom-type q x y z mux muy muz -dpd: atom-ID atom-type x y z ellipsoid: atom-ID atom-type x y z quatw quati quatj quatk full: atom-ID molecule-ID atom-type q x y z granular: atom-ID atom-type diameter density x y z diff -Naur lammps-11Nov09/src/ASPHERE/atom_vec_ellipsoid.cpp lammps-12Nov09/src/ASPHERE/atom_vec_ellipsoid.cpp --- lammps-11Nov09/src/ASPHERE/atom_vec_ellipsoid.cpp 2008-01-21 16:04:37.000000000 -0700 +++ lammps-12Nov09/src/ASPHERE/atom_vec_ellipsoid.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -35,12 +35,15 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { + molecular = 0; mass_type = 1; shape_type = 1; + comm_x_only = comm_f_only = 0; - size_comm = 7; + size_forward = 7; size_reverse = 6; size_border = 10; + size_velocity = 6; size_data_atom = 9; size_data_vel = 7; xcol_data = 3; @@ -158,6 +161,62 @@ /* ---------------------------------------------------------------------- */ +int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = quat[j][0]; + buf[m++] = quat[j][1]; + buf[m++] = quat[j][2]; + buf[m++] = quat[j][3]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = quat[j][0]; + buf[m++] = quat[j][1]; + buf[m++] = quat[j][2]; + buf[m++] = quat[j][3]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecEllipsoid::pack_comm_one(int i, double *buf) { buf[0] = quat[i][0]; @@ -188,6 +247,31 @@ /* ---------------------------------------------------------------------- */ +void AtomVecEllipsoid::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + quat[i][0] = buf[m++]; + quat[i][1] = buf[m++]; + quat[i][2] = buf[m++]; + quat[i][3] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + angmom[i][0] = buf[m++]; + angmom[i][1] = buf[m++]; + angmom[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecEllipsoid::unpack_comm_one(int i, double *buf) { quat[i][0] = buf[0]; @@ -306,6 +390,68 @@ /* ---------------------------------------------------------------------- */ +int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = quat[j][0]; + buf[m++] = quat[j][1]; + buf[m++] = quat[j][2]; + buf[m++] = quat[j][3]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = quat[j][0]; + buf[m++] = quat[j][1]; + buf[m++] = quat[j][2]; + buf[m++] = quat[j][3]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecEllipsoid::pack_border_one(int i, double *buf) { buf[0] = quat[i][0]; @@ -340,6 +486,35 @@ /* ---------------------------------------------------------------------- */ +void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + quat[i][0] = buf[m++]; + quat[i][1] = buf[m++]; + quat[i][2] = buf[m++]; + quat[i][3] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + angmom[i][0] = buf[m++]; + angmom[i][1] = buf[m++]; + angmom[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecEllipsoid::unpack_border_one(int i, double *buf) { quat[i][0] = buf[0]; diff -Naur lammps-11Nov09/src/ASPHERE/atom_vec_ellipsoid.h lammps-12Nov09/src/ASPHERE/atom_vec_ellipsoid.h --- lammps-11Nov09/src/ASPHERE/atom_vec_ellipsoid.h 2008-01-21 16:04:37.000000000 -0700 +++ lammps-12Nov09/src/ASPHERE/atom_vec_ellipsoid.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,16 +25,20 @@ void grow(int); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); int pack_comm_one(int, double *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int unpack_comm_one(int, double *); int pack_reverse(int, int, double *); int pack_reverse_one(int, double *); void unpack_reverse(int, int *, double *); int unpack_reverse_one(int, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/COLLOID/atom_vec_colloid.cpp lammps-12Nov09/src/COLLOID/atom_vec_colloid.cpp --- lammps-11Nov09/src/COLLOID/atom_vec_colloid.cpp 2009-11-06 14:02:28.000000000 -0700 +++ lammps-12Nov09/src/COLLOID/atom_vec_colloid.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -30,13 +30,16 @@ AtomVecColloid::AtomVecColloid(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { + molecular = 0; mass_type = 1; shape_type = 1; - comm_x_only = comm_f_only = 0; - ghost_velocity = 1; - size_comm = 9; + + comm_x_only = 1; + comm_f_only = 0; + size_forward = 3; size_reverse = 6; - size_border = 12; + size_border = 6; + size_velocity = 6; size_data_atom = 5; size_data_vel = 7; xcol_data = 3; @@ -117,6 +120,42 @@ buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecColloid::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; @@ -152,20 +191,23 @@ /* ---------------------------------------------------------------------- */ -int AtomVecColloid::pack_comm_one(int i, double *buf) +void AtomVecColloid::unpack_comm(int n, int first, double *buf) { - buf[0] = v[i][0]; - buf[1] = v[i][1]; - buf[2] = v[i][2]; - buf[3] = omega[i][0]; - buf[4] = omega[i][1]; - buf[5] = omega[i][2]; - return 6; + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + } } + /* ---------------------------------------------------------------------- */ -void AtomVecColloid::unpack_comm(int n, int first, double *buf) +void AtomVecColloid::unpack_comm_vel(int n, int first, double *buf) { int i,m,last; @@ -186,19 +228,6 @@ /* ---------------------------------------------------------------------- */ -int AtomVecColloid::unpack_comm_one(int i, double *buf) -{ - v[i][0] = buf[0]; - v[i][1] = buf[1]; - v[i][2] = buf[2]; - omega[i][0] = buf[3]; - omega[i][1] = buf[4]; - omega[i][2] = buf[5]; - return 6; -} - -/* ---------------------------------------------------------------------- */ - int AtomVecColloid::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -272,6 +301,48 @@ buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecColloid::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; @@ -310,20 +381,26 @@ /* ---------------------------------------------------------------------- */ -int AtomVecColloid::pack_border_one(int i, double *buf) +void AtomVecColloid::unpack_border(int n, int first, double *buf) { - buf[0] = v[i][0]; - buf[1] = v[i][1]; - buf[2] = v[i][2]; - buf[3] = omega[i][0]; - buf[4] = omega[i][1]; - buf[5] = omega[i][2]; - return 6; + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + } } /* ---------------------------------------------------------------------- */ -void AtomVecColloid::unpack_border(int n, int first, double *buf) +void AtomVecColloid::unpack_border_vel(int n, int first, double *buf) { int i,m,last; @@ -346,19 +423,6 @@ } } -/* ---------------------------------------------------------------------- */ - -int AtomVecColloid::unpack_border_one(int i, double *buf) -{ - v[i][0] = buf[0]; - v[i][1] = buf[1]; - v[i][2] = buf[2]; - omega[i][0] = buf[3]; - omega[i][1] = buf[4]; - omega[i][2] = buf[5]; - return 6; -} - /* ---------------------------------------------------------------------- pack data for atom I for sending to another proc xyz must be 1st 3 values, so comm::exchange() can test on them diff -Naur lammps-11Nov09/src/COLLOID/atom_vec_colloid.h lammps-12Nov09/src/COLLOID/atom_vec_colloid.h --- lammps-11Nov09/src/COLLOID/atom_vec_colloid.h 2009-11-06 14:02:28.000000000 -0700 +++ lammps-12Nov09/src/COLLOID/atom_vec_colloid.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,17 +25,17 @@ void grow(int); void copy(int, int); int pack_comm(int, int *, double *, int, int *); - int pack_comm_one(int, double *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); - int unpack_comm_one(int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); int pack_reverse_one(int, double *); void unpack_reverse(int, int *, double *); int unpack_reverse_one(int, double *); int pack_border(int, int *, double *, int, int *); - int pack_border_one(int, double *); + int pack_border_vel(int, int *, double *, int, int *); void unpack_border(int, int, double *); - int unpack_border_one(int, double *); + void unpack_border_vel(int, int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); int size_restart(); diff -Naur lammps-11Nov09/src/COLLOID/pair_lubricate.cpp lammps-12Nov09/src/COLLOID/pair_lubricate.cpp --- lammps-11Nov09/src/COLLOID/pair_lubricate.cpp 2009-08-08 14:06:53.000000000 -0600 +++ lammps-12Nov09/src/COLLOID/pair_lubricate.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -390,6 +390,8 @@ void PairLubricate::init_style() { + if (comm->ghost_velocity == 0) + error->all("Pair lubricate requires ghost atoms store velocity"); if (!atom->torque_flag || !atom->avec->shape_type) error->all("Pair lubricate requires atom attributes torque and shape"); if (atom->radius_flag || atom->rmass_flag) diff -Naur lammps-11Nov09/src/DIPOLE/atom_vec_dipole.cpp lammps-12Nov09/src/DIPOLE/atom_vec_dipole.cpp --- lammps-11Nov09/src/DIPOLE/atom_vec_dipole.cpp 2008-01-21 16:04:37.000000000 -0700 +++ lammps-12Nov09/src/DIPOLE/atom_vec_dipole.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -29,13 +29,16 @@ AtomVecDipole::AtomVecDipole(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { + molecular = 0; mass_type = 1; shape_type = 1; dipole_type = 1; + comm_x_only = comm_f_only = 0; - size_comm = 6; + size_forward = 6; size_reverse = 6; size_border = 10; + size_velocity = 6; size_data_atom = 9; size_data_vel = 7; xcol_data = 4; @@ -153,6 +156,60 @@ /* ---------------------------------------------------------------------- */ +int AtomVecDipole::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecDipole::pack_comm_one(int i, double *buf) { buf[0] = mu[i][0]; @@ -181,6 +238,30 @@ /* ---------------------------------------------------------------------- */ +void AtomVecDipole::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + mu[i][0] = buf[m++]; + mu[i][1] = buf[m++]; + mu[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecDipole::unpack_comm_one(int i, double *buf) { mu[i][0] = buf[0]; @@ -298,6 +379,68 @@ /* ---------------------------------------------------------------------- */ +int AtomVecDipole::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = mu[j][0]; + buf[m++] = mu[j][1]; + buf[m++] = mu[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecDipole::pack_border_one(int i, double *buf) { buf[0] = q[i]; @@ -332,6 +475,35 @@ /* ---------------------------------------------------------------------- */ +void AtomVecDipole::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + q[i] = buf[m++]; + mu[i][0] = buf[m++]; + mu[i][1] = buf[m++]; + mu[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecDipole::unpack_border_one(int i, double *buf) { q[i] = buf[0]; diff -Naur lammps-11Nov09/src/DIPOLE/atom_vec_dipole.h lammps-12Nov09/src/DIPOLE/atom_vec_dipole.h --- lammps-11Nov09/src/DIPOLE/atom_vec_dipole.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/DIPOLE/atom_vec_dipole.h 2009-11-09 11:20:20.000000000 -0700 @@ -24,16 +24,20 @@ void grow(int); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); int pack_comm_one(int, double *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int unpack_comm_one(int, double *); int pack_reverse(int, int, double *); int pack_reverse_one(int, double *); void unpack_reverse(int, int *, double *); int unpack_reverse_one(int, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/DPD/Install.csh lammps-12Nov09/src/DPD/Install.csh --- lammps-11Nov09/src/DPD/Install.csh 2007-01-29 17:22:05.000000000 -0700 +++ lammps-12Nov09/src/DPD/Install.csh 1969-12-31 17:00:00.000000000 -0700 @@ -1,24 +0,0 @@ -# Install/unInstall package classes in LAMMPS - -if ($1 == 1) then - - cp style_dpd.h .. - - cp atom_vec_dpd.cpp .. - cp pair_dpd.cpp .. - - cp atom_vec_dpd.h .. - cp pair_dpd.h .. - -else if ($1 == 0) then - - rm ../style_dpd.h - touch ../style_dpd.h - - rm ../atom_vec_dpd.cpp - rm ../pair_dpd.cpp - - rm ../atom_vec_dpd.h - rm ../pair_dpd.h - -endif diff -Naur lammps-11Nov09/src/DPD/atom_vec_dpd.cpp lammps-12Nov09/src/DPD/atom_vec_dpd.cpp --- lammps-11Nov09/src/DPD/atom_vec_dpd.cpp 2008-02-08 12:47:58.000000000 -0700 +++ lammps-12Nov09/src/DPD/atom_vec_dpd.cpp 1969-12-31 17:00:00.000000000 -0700 @@ -1,211 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "stdlib.h" -#include "atom_vec_dpd.h" -#include "atom.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define DELTA 10000 - -/* ---------------------------------------------------------------------- */ - -AtomVecDPD::AtomVecDPD(LAMMPS *lmp, int narg, char **arg) : - AtomVecAtomic(lmp, narg, arg) -{ - mass_type = 1; - comm_x_only = 0; - ghost_velocity = 1; - size_comm = 6; - size_reverse = 3; - size_border = 9; - size_data_atom = 5; - size_data_vel = 4; - xcol_data = 3; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecDPD::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; - dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; - dz = pbc[2]*domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecDPD::pack_comm_one(int i, double *buf) -{ - buf[0] = v[i][0]; - buf[1] = v[i][1]; - buf[2] = v[i][2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecDPD::unpack_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecDPD::unpack_comm_one(int i, double *buf) -{ - v[i][0] = buf[0]; - v[i][1] = buf[1]; - v[i][2] = buf[2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecDPD::pack_border(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = tag[j]; - buf[m++] = type[j]; - buf[m++] = mask[j]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = tag[j]; - buf[m++] = type[j]; - buf[m++] = mask[j]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecDPD::pack_border_one(int i, double *buf) -{ - buf[0] = v[i][0]; - buf[1] = v[i][1]; - buf[2] = v[i][2]; - return 3; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecDPD::unpack_border(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = static_cast (buf[m++]); - type[i] = static_cast (buf[m++]); - mask[i] = static_cast (buf[m++]); - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecDPD::unpack_border_one(int i, double *buf) -{ - v[i][0] = buf[0]; - v[i][1] = buf[1]; - v[i][2] = buf[2]; - return 3; -} diff -Naur lammps-11Nov09/src/DPD/atom_vec_dpd.h lammps-12Nov09/src/DPD/atom_vec_dpd.h --- lammps-11Nov09/src/DPD/atom_vec_dpd.h 2007-03-15 15:49:48.000000000 -0600 +++ lammps-12Nov09/src/DPD/atom_vec_dpd.h 1969-12-31 17:00:00.000000000 -0700 @@ -1,37 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 ATOM_VEC_DPD_H -#define ATOM_VEC_DPD_H - -#include "atom_vec_atomic.h" - -namespace LAMMPS_NS { - -class AtomVecDPD : public AtomVecAtomic { - public: - AtomVecDPD(class LAMMPS *, int, char **); - void zero_ghost(int, int); - int pack_comm(int, int *, double *, int, int *); - int pack_comm_one(int, double *); - void unpack_comm(int, int, double *); - int unpack_comm_one(int, double *); - int pack_border(int, int *, double *, int, int *); - int pack_border_one(int, double *); - void unpack_border(int, int, double *); - int unpack_border_one(int, double *); -}; - -} - -#endif diff -Naur lammps-11Nov09/src/DPD/pair_dpd.cpp lammps-12Nov09/src/DPD/pair_dpd.cpp --- lammps-11Nov09/src/DPD/pair_dpd.cpp 2009-11-02 11:25:41.000000000 -0700 +++ lammps-12Nov09/src/DPD/pair_dpd.cpp 1969-12-31 17:00:00.000000000 -0700 @@ -1,388 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Kurt Smith (U Pittsburgh) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "pair_dpd.h" -#include "atom.h" -#include "atom_vec.h" -#include "comm.h" -#include "update.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "random_mars.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - -#define EPSILON 1.0e-10 - -/* ---------------------------------------------------------------------- */ - -PairDPD::PairDPD(LAMMPS *lmp) : Pair(lmp) -{ - random = NULL; -} - -/* ---------------------------------------------------------------------- */ - -PairDPD::~PairDPD() -{ - if (allocated) { - memory->destroy_2d_int_array(setflag); - memory->destroy_2d_double_array(cutsq); - - memory->destroy_2d_double_array(cut); - memory->destroy_2d_double_array(a0); - memory->destroy_2d_double_array(gamma); - memory->destroy_2d_double_array(sigma); - } - - if (random) delete random; -} - -/* ---------------------------------------------------------------------- */ - -void PairDPD::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double vxtmp,vytmp,vztmp,delvx,delvy,delvz; - double rsq,r,rinv,dot,wd,randnum,factor_dpd; - int *ilist,*jlist,*numneigh,**firstneigh; - - evdwl = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - int *type = atom->type; - int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double dtinvsqrt = 1.0/sqrt(update->dt); - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - vxtmp = v[i][0]; - vytmp = v[i][1]; - vztmp = v[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - - if (j < nall) factor_dpd = 1.0; - else { - factor_dpd = special_lj[j/nall]; - j %= nall; - } - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r = sqrt(rsq); - if (r < EPSILON) continue; // r can be 0.0 in DPD systems - rinv = 1.0/r; - delvx = vxtmp - v[j][0]; - delvy = vytmp - v[j][1]; - delvz = vztmp - v[j][2]; - dot = delx*delvx + dely*delvy + delz*delvz; - wd = 1.0 - r/cut[itype][jtype]; - randnum = random->gaussian(); - - // conservative force = a0 * wd - // drag force = -gamma * wd^2 * (delx dot delv) / r - // random force = sigma * wd * rnd * dtinvsqrt; - - fpair = a0[itype][jtype]*wd; - fpair -= gamma[itype][jtype]*wd*wd*dot*rinv; - fpair += sigma[itype][jtype]*wd*randnum*dtinvsqrt; - fpair *= factor_dpd*rinv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; - } - - if (eflag) { - evdwl = -a0[itype][jtype] * r * (1.0 - 0.5*r/cut[itype][jtype]); - evdwl *= factor_dpd; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); - } - } - } - - if (vflag_fdotr) virial_compute(); -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairDPD::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); - - cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); - a0 = memory->create_2d_double_array(n+1,n+1,"pair:a0"); - gamma = memory->create_2d_double_array(n+1,n+1,"pair:gamma"); - sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairDPD::settings(int narg, char **arg) -{ - if (narg != 3) error->all("Illegal pair_style command"); - - temperature = force->numeric(arg[0]); - cut_global = force->numeric(arg[1]); - seed = force->inumeric(arg[2]); - - // initialize Marsaglia RNG with processor-unique seed - - if (seed <= 0) error->all("Illegal pair_style command"); - delete random; - random = new RanMars(lmp,seed + comm->me); - - // reset cutoffs that have been explicitly set - - if (allocated) { - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i+1; j <= atom->ntypes; j++) - if (setflag[i][j]) cut[i][j] = cut_global; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairDPD::coeff(int narg, char **arg) -{ - if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients"); - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(arg[0],atom->ntypes,ilo,ihi); - force->bounds(arg[1],atom->ntypes,jlo,jhi); - - double a0_one = force->numeric(arg[2]); - double gamma_one = force->numeric(arg[3]); - - double cut_one = cut_global; - if (narg == 5) cut_one = force->numeric(arg[4]); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - a0[i][j] = a0_one; - gamma[i][j] = gamma_one; - cut[i][j] = cut_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all("Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairDPD::init_style() -{ - if (atom->avec->ghost_velocity == 0) - error->all("Pair dpd requires ghost atoms store velocity"); - - // if newton off, forces between atoms ij will be double computed - // using different random numbers - - if (force->newton_pair == 0 && comm->me == 0) error->warning( - "Pair dpd needs newton pair on for momentum conservation"); - - int irequest = neighbor->request(this); -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairDPD::init_one(int i, int j) -{ - if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); - - sigma[i][j] = sqrt(2.0*force->boltz*temperature*gamma[i][j]); - - cut[j][i] = cut[i][j]; - a0[j][i] = a0[i][j]; - gamma[j][i] = gamma[i][j]; - sigma[j][i] = sigma[i][j]; - - return cut[i][j]; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairDPD::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&a0[i][j],sizeof(double),1,fp); - fwrite(&gamma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairDPD::read_restart(FILE *fp) -{ - read_restart_settings(fp); - - allocate(); - - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - fread(&a0[i][j],sizeof(double),1,fp); - fread(&gamma[i][j],sizeof(double),1,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&gamma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairDPD::write_restart_settings(FILE *fp) -{ - fwrite(&temperature,sizeof(double),1,fp); - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&seed,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairDPD::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&temperature,sizeof(double),1,fp); - fread(&cut_global,sizeof(double),1,fp); - fread(&seed,sizeof(int),1,fp); - fread(&mix_flag,sizeof(int),1,fp); - } - MPI_Bcast(&temperature,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&seed,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - - // initialize Marsaglia RNG with processor-unique seed - // same seed that pair_style command initially specified - - if (random) delete random; - random = new RanMars(lmp,seed + comm->me); -} - -/* ---------------------------------------------------------------------- */ - -double PairDPD::single(int i, int j, int itype, int jtype, double rsq, - double factor_coul, double factor_dpd, double &fforce) -{ - double r,rinv,wd,phi; - - r = sqrt(rsq); - if (r < EPSILON) { - fforce = 0.0; - return 0.0; - } - - rinv = 1.0/r; - wd = 1.0 - r/cut[itype][jtype]; - fforce = a0[itype][jtype]*wd * factor_dpd*rinv; - - phi = -a0[itype][jtype] * r * (1.0 - 0.5*r/cut[itype][jtype]); - return factor_dpd*phi; -} diff -Naur lammps-11Nov09/src/DPD/pair_dpd.h lammps-12Nov09/src/DPD/pair_dpd.h --- lammps-11Nov09/src/DPD/pair_dpd.h 2008-01-02 12:24:46.000000000 -0700 +++ lammps-12Nov09/src/DPD/pair_dpd.h 1969-12-31 17:00:00.000000000 -0700 @@ -1,49 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 PAIR_DPD_H -#define PAIR_DPD_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairDPD : public Pair { - public: - PairDPD(class LAMMPS *); - ~PairDPD(); - void compute(int, int); - void settings(int, char **); - void coeff(int, char **); - void init_style(); - double init_one(int, int); - void write_restart(FILE *); - void read_restart(FILE *); - void write_restart_settings(FILE *); - void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); - - private: - double cut_global,temperature; - int seed; - double **cut; - double **a0,**gamma; - double **sigma; - class RanMars *random; - - void allocate(); -}; - -} - -#endif diff -Naur lammps-11Nov09/src/DPD/style_dpd.h lammps-12Nov09/src/DPD/style_dpd.h --- lammps-11Nov09/src/DPD/style_dpd.h 2007-01-29 17:22:05.000000000 -0700 +++ lammps-12Nov09/src/DPD/style_dpd.h 1969-12-31 17:00:00.000000000 -0700 @@ -1,28 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 AtomInclude -#include "atom_vec_dpd.h" -#endif - -#ifdef AtomClass -AtomStyle(dpd,AtomVecDPD) -#endif - -#ifdef PairInclude -#include "pair_dpd.h" -#endif - -#ifdef PairClass -PairStyle(dpd,PairDPD) -#endif diff -Naur lammps-11Nov09/src/GRANULAR/atom_vec_granular.cpp lammps-12Nov09/src/GRANULAR/atom_vec_granular.cpp --- lammps-11Nov09/src/GRANULAR/atom_vec_granular.cpp 2009-06-26 12:23:16.000000000 -0600 +++ lammps-12Nov09/src/GRANULAR/atom_vec_granular.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -31,11 +31,14 @@ AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { - comm_x_only = comm_f_only = 0; - ghost_velocity = 1; - size_comm = 9; + molecular = 0; + + comm_x_only = 1; + comm_f_only = 0; + size_forward = 3; size_reverse = 6; - size_border = 14; + size_border = 8; + size_velocity = 6; size_data_atom = 7; size_data_vel = 7; xcol_data = 5; @@ -129,6 +132,42 @@ buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; @@ -164,20 +203,22 @@ /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_comm_one(int i, double *buf) +void AtomVecGranular::unpack_comm(int n, int first, double *buf) { - buf[0] = v[i][0]; - buf[1] = v[i][1]; - buf[2] = v[i][2]; - buf[3] = omega[i][0]; - buf[4] = omega[i][1]; - buf[5] = omega[i][2]; - return 6; + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + } } /* ---------------------------------------------------------------------- */ -void AtomVecGranular::unpack_comm(int n, int first, double *buf) +void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf) { int i,m,last; @@ -198,19 +239,6 @@ /* ---------------------------------------------------------------------- */ -int AtomVecGranular::unpack_comm_one(int i, double *buf) -{ - v[i][0] = buf[0]; - v[i][1] = buf[1]; - v[i][2] = buf[2]; - omega[i][0] = buf[3]; - omega[i][1] = buf[4]; - omega[i][2] = buf[5]; - return 6; -} - -/* ---------------------------------------------------------------------- */ - int AtomVecGranular::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -284,11 +312,57 @@ buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; + buf[m++] = radius[j]; + buf[m++] = rmass[j]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = radius[j]; + buf[m++] = rmass[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGranular::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = radius[j]; + buf[m++] = rmass[j]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; - buf[m++] = radius[j]; - buf[m++] = rmass[j]; buf[m++] = omega[j][0]; buf[m++] = omega[j][1]; buf[m++] = omega[j][2]; @@ -311,11 +385,11 @@ buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; + buf[m++] = radius[j]; + buf[m++] = rmass[j]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; - buf[m++] = radius[j]; - buf[m++] = rmass[j]; buf[m++] = omega[j][0]; buf[m++] = omega[j][1]; buf[m++] = omega[j][2]; @@ -328,15 +402,9 @@ int AtomVecGranular::pack_border_one(int i, double *buf) { - buf[0] = v[i][0]; - buf[1] = v[i][1]; - buf[2] = v[i][2]; - buf[3] = radius[i]; - buf[4] = rmass[i]; - buf[5] = omega[i][0]; - buf[6] = omega[i][1]; - buf[7] = omega[i][2]; - return 8; + buf[0] = radius[i]; + buf[1] = rmass[i]; + return 2; } /* ---------------------------------------------------------------------- */ @@ -355,11 +423,33 @@ tag[i] = static_cast (buf[m++]); type[i] = static_cast (buf[m++]); mask[i] = static_cast (buf[m++]); + radius[i] = buf[m++]; + rmass[i] = buf[m++]; + } +} + + +/* ---------------------------------------------------------------------- */ + +void AtomVecGranular::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + radius[i] = buf[m++]; + rmass[i] = buf[m++]; v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; - radius[i] = buf[m++]; - rmass[i] = buf[m++]; omega[i][0] = buf[m++]; omega[i][1] = buf[m++]; omega[i][2] = buf[m++]; @@ -370,15 +460,9 @@ int AtomVecGranular::unpack_border_one(int i, double *buf) { - v[i][0] = buf[0]; - v[i][1] = buf[1]; - v[i][2] = buf[2]; - radius[i] = buf[3]; - rmass[i] = buf[4]; - omega[i][0] = buf[5]; - omega[i][1] = buf[6]; - omega[i][2] = buf[7]; - return 8; + radius[i] = buf[0]; + rmass[i] = buf[1]; + return 2; } /* ---------------------------------------------------------------------- diff -Naur lammps-11Nov09/src/GRANULAR/atom_vec_granular.h lammps-12Nov09/src/GRANULAR/atom_vec_granular.h --- lammps-11Nov09/src/GRANULAR/atom_vec_granular.h 2008-03-17 17:37:19.000000000 -0600 +++ lammps-12Nov09/src/GRANULAR/atom_vec_granular.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,16 +25,18 @@ void grow(int); void copy(int, int); int pack_comm(int, int *, double *, int, int *); - int pack_comm_one(int, double *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); - int unpack_comm_one(int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); int pack_reverse_one(int, double *); void unpack_reverse(int, int *, double *); int unpack_reverse_one(int, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/GRANULAR/pair_gran_hooke_history.cpp lammps-12Nov09/src/GRANULAR/pair_gran_hooke_history.cpp --- lammps-11Nov09/src/GRANULAR/pair_gran_hooke_history.cpp 2009-08-08 14:06:53.000000000 -0600 +++ lammps-12Nov09/src/GRANULAR/pair_gran_hooke_history.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -343,7 +343,7 @@ if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag) error->all("Pair granular requires atom attributes radius, omega, torque"); - if (atom->avec->ghost_velocity == 0) + if (comm->ghost_velocity == 0) error->all("Pair granular requires ghost atoms store velocity"); // need a half neigh list and optionally a granular history neigh list diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_angle.cpp lammps-12Nov09/src/MOLECULE/atom_vec_angle.cpp --- lammps-11Nov09/src/MOLECULE/atom_vec_angle.cpp 2009-07-21 09:00:56.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_angle.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -33,12 +33,14 @@ AtomVec(lmp, narg, arg) { molecular = 1; - bonds_allow = 1; - angles_allow = 1; + bonds_allow = angles_allow = 1; mass_type = 1; - size_comm = 3; + + comm_x_only = comm_f_only = 1; + size_forward = 3; size_reverse = 3; size_border = 7; + size_velocity = 3; size_data_atom = 6; size_data_vel = 4; xcol_data = 4; @@ -173,6 +175,48 @@ buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecAngle::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; } } else { if (domain->triclinic == 0) { @@ -211,6 +255,24 @@ /* ---------------------------------------------------------------------- */ +void AtomVecAngle::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecAngle::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -286,6 +348,56 @@ /* ---------------------------------------------------------------------- */ +int AtomVecAngle::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecAngle::pack_border_one(int i, double *buf) { buf[0] = molecule[i]; @@ -314,6 +426,29 @@ /* ---------------------------------------------------------------------- */ +void AtomVecAngle::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + molecule[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecAngle::unpack_border_one(int i, double *buf) { molecule[i] = static_cast (buf[0]); diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_angle.h lammps-12Nov09/src/MOLECULE/atom_vec_angle.h --- lammps-11Nov09/src/MOLECULE/atom_vec_angle.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_angle.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,12 +25,16 @@ void reset_special(); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_bond.cpp lammps-12Nov09/src/MOLECULE/atom_vec_bond.cpp --- lammps-11Nov09/src/MOLECULE/atom_vec_bond.cpp 2009-07-21 09:00:56.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_bond.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -35,9 +35,12 @@ molecular = 1; bonds_allow = 1; mass_type = 1; - size_comm = 3; + + comm_x_only = comm_f_only = 1; + size_forward = 3; size_reverse = 3; size_border = 7; + size_velocity = 3; size_data_atom = 6; size_data_vel = 4; xcol_data = 4; @@ -172,6 +175,48 @@ /* ---------------------------------------------------------------------- */ +int AtomVecBond::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecBond::unpack_comm(int n, int first, double *buf) { int i,m,last; @@ -187,6 +232,24 @@ /* ---------------------------------------------------------------------- */ +void AtomVecBond::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecBond::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -262,6 +325,56 @@ /* ---------------------------------------------------------------------- */ +int AtomVecBond::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecBond::pack_border_one(int i, double *buf) { buf[0] = molecule[i]; @@ -290,6 +403,29 @@ /* ---------------------------------------------------------------------- */ +void AtomVecBond::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + molecule[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecBond::unpack_border_one(int i, double *buf) { molecule[i] = static_cast (buf[0]); diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_bond.h lammps-12Nov09/src/MOLECULE/atom_vec_bond.h --- lammps-11Nov09/src/MOLECULE/atom_vec_bond.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_bond.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,12 +25,16 @@ void reset_special(); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_full.cpp lammps-12Nov09/src/MOLECULE/atom_vec_full.cpp --- lammps-11Nov09/src/MOLECULE/atom_vec_full.cpp 2009-07-21 09:00:56.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_full.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -33,14 +33,14 @@ AtomVec(lmp, narg, arg) { molecular = 1; - bonds_allow = 1; - angles_allow = 1; - dihedrals_allow = 1; - impropers_allow = 1; + bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 1; mass_type = 1; - size_comm = 3; + + comm_x_only = comm_f_only = 1; + size_forward = 3; size_reverse = 3; size_border = 8; + size_velocity = 3; size_data_atom = 7; size_data_vel = 4; xcol_data = 5; @@ -255,6 +255,48 @@ /* ---------------------------------------------------------------------- */ +int AtomVecFull::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecFull::unpack_comm(int n, int first, double *buf) { int i,m,last; @@ -270,6 +312,24 @@ /* ---------------------------------------------------------------------- */ +void AtomVecFull::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecFull::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -347,6 +407,58 @@ /* ---------------------------------------------------------------------- */ +int AtomVecFull::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecFull::pack_border_one(int i, double *buf) { buf[0] = q[i]; @@ -377,6 +489,30 @@ /* ---------------------------------------------------------------------- */ +void AtomVecFull::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + q[i] = buf[m++]; + molecule[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecFull::unpack_border_one(int i, double *buf) { q[i] = buf[0]; diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_full.h lammps-12Nov09/src/MOLECULE/atom_vec_full.h --- lammps-11Nov09/src/MOLECULE/atom_vec_full.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_full.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,12 +25,16 @@ void reset_special(); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_molecular.cpp lammps-12Nov09/src/MOLECULE/atom_vec_molecular.cpp --- lammps-11Nov09/src/MOLECULE/atom_vec_molecular.cpp 2009-07-21 09:00:56.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_molecular.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -33,14 +33,14 @@ AtomVec(lmp, narg, arg) { molecular = 1; - bonds_allow = 1; - angles_allow = 1; - dihedrals_allow = 1; - impropers_allow = 1; + bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 1; mass_type = 1; - size_comm = 3; + + comm_x_only = comm_f_only = 1; + size_forward = 3; size_reverse = 3; size_border = 7; + size_velocity = 3; size_data_atom = 6; size_data_vel = 4; xcol_data = 4; @@ -252,6 +252,48 @@ /* ---------------------------------------------------------------------- */ +int AtomVecMolecular::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecMolecular::unpack_comm(int n, int first, double *buf) { int i,m,last; @@ -267,6 +309,24 @@ /* ---------------------------------------------------------------------- */ +void AtomVecMolecular::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecMolecular::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -342,6 +402,56 @@ /* ---------------------------------------------------------------------- */ +int AtomVecMolecular::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = molecule[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecMolecular::pack_border_one(int i, double *buf) { buf[0] = molecule[i]; @@ -370,6 +480,29 @@ /* ---------------------------------------------------------------------- */ +void AtomVecMolecular::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + molecule[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecMolecular::unpack_border_one(int i, double *buf) { molecule[i] = static_cast (buf[0]); diff -Naur lammps-11Nov09/src/MOLECULE/atom_vec_molecular.h lammps-12Nov09/src/MOLECULE/atom_vec_molecular.h --- lammps-11Nov09/src/MOLECULE/atom_vec_molecular.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/MOLECULE/atom_vec_molecular.h 2009-11-09 11:20:20.000000000 -0700 @@ -25,12 +25,16 @@ void reset_special(); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/Makefile lammps-12Nov09/src/Makefile --- lammps-11Nov09/src/Makefile 2009-11-06 09:48:55.000000000 -0700 +++ lammps-12Nov09/src/Makefile 2009-11-09 11:20:20.000000000 -0700 @@ -13,7 +13,7 @@ # Package variables -PACKAGE = asphere class2 colloid dipole dpd dsmc gpu granular \ +PACKAGE = asphere class2 colloid dipole dsmc gpu granular \ kspace manybody meam molecule opt peri poems prd reax xtc PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \ diff -Naur lammps-11Nov09/src/PERI/atom_vec_peri.cpp lammps-12Nov09/src/PERI/atom_vec_peri.cpp --- lammps-11Nov09/src/PERI/atom_vec_peri.cpp 2009-05-19 09:02:32.000000000 -0600 +++ lammps-12Nov09/src/PERI/atom_vec_peri.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -32,12 +32,16 @@ /* ---------------------------------------------------------------------- */ AtomVecPeri::AtomVecPeri(LAMMPS *lmp, int narg, char **arg) : -AtomVec(lmp, narg, arg) + AtomVec(lmp, narg, arg) { + molecular = 0; + comm_x_only = 0; - size_comm = 4; + comm_f_only = 1; + size_forward = 4; size_reverse = 3; size_border = 11; + size_velocity = 3; size_data_atom = 7; size_data_vel = 4; xcol_data = 5; @@ -153,6 +157,51 @@ /* ---------------------------------------------------------------------- */ +int AtomVecPeri::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) + +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = s0[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = s0[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecPeri::pack_comm_one(int i, double *buf) { buf[0] = s0[i]; @@ -177,6 +226,25 @@ /* ---------------------------------------------------------------------- */ +void AtomVecPeri::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + s0[i] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecPeri::unpack_comm_one(int i, double *buf) { s0[i] = buf[0]; @@ -268,6 +336,64 @@ /* ---------------------------------------------------------------------- */ +int AtomVecPeri::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = vfrac[j]; + buf[m++] = s0[j]; + buf[m++] = x0[j][0]; + buf[m++] = x0[j][1]; + buf[m++] = x0[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = vfrac[j]; + buf[m++] = s0[j]; + buf[m++] = x0[j][0]; + buf[m++] = x0[j][1]; + buf[m++] = x0[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecPeri::pack_border_one(int i, double *buf) { buf[0] = vfrac[i]; @@ -303,6 +429,33 @@ } /* ---------------------------------------------------------------------- */ + +void AtomVecPeri::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + vfrac[i] = buf[m++]; + s0[i] = buf[m++]; + x0[i][0] = buf[m++]; + x0[i][1] = buf[m++]; + x0[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ int AtomVecPeri::unpack_border_one(int i, double *buf) { diff -Naur lammps-11Nov09/src/PERI/atom_vec_peri.h lammps-12Nov09/src/PERI/atom_vec_peri.h --- lammps-11Nov09/src/PERI/atom_vec_peri.h 2008-07-25 08:41:48.000000000 -0600 +++ lammps-12Nov09/src/PERI/atom_vec_peri.h 2009-11-09 11:20:20.000000000 -0700 @@ -24,14 +24,18 @@ void grow(int); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); int pack_comm_one(int, double *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int unpack_comm_one(int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/atom_vec.cpp lammps-12Nov09/src/atom_vec.cpp --- lammps-11Nov09/src/atom_vec.cpp 2008-02-08 12:47:58.000000000 -0700 +++ lammps-12Nov09/src/atom_vec.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -22,12 +22,8 @@ AtomVec::AtomVec(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) { nmax = 0; - molecular = 0; bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 0; mass_type = shape_type = dipole_type = 0; - comm_x_only = comm_f_only = 1; - ghost_velocity = 0; - size_comm = size_reverse = size_border = 0; } /* ---------------------------------------------------------------------- diff -Naur lammps-11Nov09/src/atom_vec.h lammps-12Nov09/src/atom_vec.h --- lammps-11Nov09/src/atom_vec.h 2008-02-08 12:47:58.000000000 -0700 +++ lammps-12Nov09/src/atom_vec.h 2009-11-09 11:20:20.000000000 -0700 @@ -26,13 +26,14 @@ int mass_type; // 1 if per-type masses int shape_type; // 1 if per-type shape array int dipole_type; // 1 if per-type dipole moments + int comm_x_only; // 1 if only exchange x in forward comm int comm_f_only; // 1 if only exchange f in reverse comm - int ghost_velocity; // 1 if ghost atoms store velocity - int size_comm; // # of values per atom in comm + int size_forward; // # of values per atom in comm int size_reverse; // # in reverse comm int size_border; // # in border comm + int size_velocity; // # of velocity based quantities int size_data_atom; // number of values in Atom line int size_data_vel; // number of values in Velocity line int xcol_data; // column (1-N) where x is in Atom line @@ -46,8 +47,10 @@ virtual void copy(int, int) = 0; virtual int pack_comm(int, int *, double *, int, int *) = 0; + virtual int pack_comm_vel(int, int *, double *, int, int *) = 0; virtual int pack_comm_one(int, double *) {return 0;} virtual void unpack_comm(int, int, double *) = 0; + virtual void unpack_comm_vel(int, int, double *) = 0; virtual int unpack_comm_one(int, double *) {return 0;} virtual int pack_reverse(int, int, double *) = 0; @@ -56,8 +59,10 @@ virtual int unpack_reverse_one(int, double *) {return 0;} virtual int pack_border(int, int *, double *, int, int *) = 0; + virtual int pack_border_vel(int, int *, double *, int, int *) = 0; virtual int pack_border_one(int, double *) {return 0;} virtual void unpack_border(int, int, double *) = 0; + virtual void unpack_border_vel(int, int, double *) = 0; virtual int unpack_border_one(int, double *) {return 0;} virtual int pack_exchange(int, double *) = 0; diff -Naur lammps-11Nov09/src/atom_vec_atomic.cpp lammps-12Nov09/src/atom_vec_atomic.cpp --- lammps-11Nov09/src/atom_vec_atomic.cpp 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/atom_vec_atomic.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -29,10 +29,14 @@ AtomVecAtomic::AtomVecAtomic(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { + molecular = 0; mass_type = 1; - size_comm = 3; + + comm_x_only = comm_f_only = 1; + size_forward = 3; size_reverse = 3; size_border = 6; + size_velocity = 3; size_data_atom = 5; size_data_vel = 4; xcol_data = 3; @@ -125,6 +129,48 @@ /* ---------------------------------------------------------------------- */ +int AtomVecAtomic::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecAtomic::unpack_comm(int n, int first, double *buf) { int i,m,last; @@ -140,6 +186,24 @@ /* ---------------------------------------------------------------------- */ +void AtomVecAtomic::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecAtomic::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -213,6 +277,54 @@ /* ---------------------------------------------------------------------- */ +int AtomVecAtomic::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecAtomic::unpack_border(int n, int first, double *buf) { int i,m,last; @@ -230,6 +342,28 @@ } } +/* ---------------------------------------------------------------------- */ + +void AtomVecAtomic::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + /* ---------------------------------------------------------------------- pack data for atom I for sending to another proc xyz must be 1st 3 values, so comm::exchange() can test on them diff -Naur lammps-11Nov09/src/atom_vec_atomic.h lammps-12Nov09/src/atom_vec_atomic.h --- lammps-11Nov09/src/atom_vec_atomic.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/atom_vec_atomic.h 2009-11-09 11:20:20.000000000 -0700 @@ -21,15 +21,19 @@ class AtomVecAtomic : public AtomVec { public: AtomVecAtomic(class LAMMPS *, int, char **); - virtual ~AtomVecAtomic() {} + ~AtomVecAtomic() {} void grow(int); void copy(int, int); - virtual int pack_comm(int, int *, double *, int, int *); - virtual void unpack_comm(int, int, double *); + int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); - virtual int pack_border(int, int *, double *, int, int *); - virtual void unpack_border(int, int, double *); + int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); + void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); int size_restart(); @@ -40,7 +44,7 @@ int data_atom_hybrid(int, char **); double memory_usage(); - protected: + private: int *tag,*type,*mask,*image; double **x,**v,**f; }; diff -Naur lammps-11Nov09/src/atom_vec_charge.cpp lammps-12Nov09/src/atom_vec_charge.cpp --- lammps-11Nov09/src/atom_vec_charge.cpp 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/atom_vec_charge.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -29,10 +29,14 @@ AtomVecCharge::AtomVecCharge(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { + molecular = 0; mass_type = 1; - size_comm = 3; + + comm_x_only = comm_f_only = 1; + size_forward = 3; size_reverse = 3; size_border = 7; + size_velocity = 3; size_data_atom = 6; size_data_vel = 4; xcol_data = 4; @@ -132,6 +136,48 @@ /* ---------------------------------------------------------------------- */ +int AtomVecCharge::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecCharge::unpack_comm(int n, int first, double *buf) { int i,m,last; @@ -147,6 +193,24 @@ /* ---------------------------------------------------------------------- */ +void AtomVecCharge::unpack_comm_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecCharge::pack_reverse(int n, int first, double *buf) { int i,m,last; @@ -222,6 +286,56 @@ /* ---------------------------------------------------------------------- */ +int AtomVecCharge::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + int AtomVecCharge::pack_border_one(int i, double *buf) { buf[0] = q[i]; @@ -250,6 +364,29 @@ /* ---------------------------------------------------------------------- */ +void AtomVecCharge::unpack_border_vel(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + q[i] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecCharge::unpack_border_one(int i, double *buf) { q[i] = buf[0]; diff -Naur lammps-11Nov09/src/atom_vec_charge.h lammps-12Nov09/src/atom_vec_charge.h --- lammps-11Nov09/src/atom_vec_charge.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/atom_vec_charge.h 2009-11-09 11:20:20.000000000 -0700 @@ -24,12 +24,16 @@ void grow(int); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); int pack_border_one(int, double *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int unpack_border_one(int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); diff -Naur lammps-11Nov09/src/atom_vec_hybrid.cpp lammps-12Nov09/src/atom_vec_hybrid.cpp --- lammps-11Nov09/src/atom_vec_hybrid.cpp 2008-10-24 10:04:25.000000000 -0600 +++ lammps-12Nov09/src/atom_vec_hybrid.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -57,7 +57,7 @@ // hybrid settings are MAX or MIN of sub-style settings // hybrid sizes are minimial values plus extra values for each sub-style - size_comm = 3; + size_forward = 3; size_reverse = 3; size_border = 6; size_data_atom = 5; @@ -73,16 +73,19 @@ mass_type = MAX(mass_type,styles[k]->mass_type); shape_type = MAX(shape_type,styles[k]->shape_type); dipole_type = MAX(dipole_type,styles[k]->dipole_type); + comm_x_only = MIN(comm_x_only,styles[k]->comm_x_only); comm_f_only = MIN(comm_f_only,styles[k]->comm_f_only); - ghost_velocity = MAX(ghost_velocity,styles[k]->ghost_velocity); - - size_comm += styles[k]->size_comm - 3; + size_forward += styles[k]->size_forward - 3; size_reverse += styles[k]->size_reverse - 3; size_border += styles[k]->size_border - 6; size_data_atom += styles[k]->size_data_atom - 5; size_data_vel += styles[k]->size_data_vel - 4; } + + size_velocity = 3; + if (atom->omega_flag) size_velocity += 3; + if (atom->angmom_flag) size_velocity += 3; } /* ---------------------------------------------------------------------- */ @@ -122,6 +125,9 @@ v = atom->v; f = atom->f; + omega = atom->omega; + angmom = atom->angmom; + if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); @@ -192,6 +198,74 @@ /* ---------------------------------------------------------------------- */ +int AtomVecHybrid::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,k,m; + double dx,dy,dz; + int omega_flag = atom->omega_flag; + int angmom_flag = atom->angmom_flag; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + if (omega_flag) { + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + if (angmom_flag) { + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + for (k = 0; k < nstyles; k++) + m += styles[k]->pack_comm_one(j,&buf[m]); + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + if (omega_flag) { + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + if (angmom_flag) { + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + for (k = 0; k < nstyles; k++) + m += styles[k]->pack_comm_one(j,&buf[m]); + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecHybrid::unpack_comm(int n, int first, double *buf) { int i,k,last; @@ -209,6 +283,38 @@ /* ---------------------------------------------------------------------- */ +void AtomVecHybrid::unpack_comm_vel(int n, int first, double *buf) +{ + int i,k,last; + int omega_flag = atom->omega_flag; + int angmom_flag = atom->angmom_flag; + + int m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + if (omega_flag) { + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; + } + if (angmom_flag) { + angmom[i][0] = buf[m++]; + angmom[i][1] = buf[m++]; + angmom[i][2] = buf[m++]; + } + for (k = 0; k < nstyles; k++) + m += styles[k]->unpack_comm_one(i,&buf[m]); + } +} + +/* ---------------------------------------------------------------------- */ + int AtomVecHybrid::pack_reverse(int n, int first, double *buf) { int i,k,last; @@ -290,6 +396,80 @@ /* ---------------------------------------------------------------------- */ +int AtomVecHybrid::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,k,m; + double dx,dy,dz; + int omega_flag = atom->omega_flag; + int angmom_flag = atom->angmom_flag; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + if (omega_flag) { + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + if (angmom_flag) { + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + for (k = 0; k < nstyles; k++) + m += styles[k]->pack_border_one(j,&buf[m]); + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + if (omega_flag) { + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + if (angmom_flag) { + buf[m++] = angmom[j][0]; + buf[m++] = angmom[j][1]; + buf[m++] = angmom[j][2]; + } + for (k = 0; k < nstyles; k++) + m += styles[k]->pack_border_one(j,&buf[m]); + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + void AtomVecHybrid::unpack_border(int n, int first, double *buf) { int i,k,m,last; @@ -309,6 +489,42 @@ } } +/* ---------------------------------------------------------------------- */ + +void AtomVecHybrid::unpack_border_vel(int n, int first, double *buf) +{ + int i,k,m,last; + int omega_flag = atom->omega_flag; + int angmom_flag = atom->angmom_flag; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + if (omega_flag) { + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; + } + if (angmom_flag) { + angmom[i][0] = buf[m++]; + angmom[i][1] = buf[m++]; + angmom[i][2] = buf[m++]; + } + for (k = 0; k < nstyles; k++) + m += styles[k]->unpack_border_one(i,&buf[m]); + } +} + /* ---------------------------------------------------------------------- pack data for atom I for sending to another proc pack each sub-style one after the other diff -Naur lammps-11Nov09/src/atom_vec_hybrid.h lammps-12Nov09/src/atom_vec_hybrid.h --- lammps-11Nov09/src/atom_vec_hybrid.h 2007-10-04 11:57:04.000000000 -0600 +++ lammps-12Nov09/src/atom_vec_hybrid.h 2009-11-09 11:20:20.000000000 -0700 @@ -30,11 +30,15 @@ void reset_special(); void copy(int, int); int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); int size_restart(); @@ -49,6 +53,7 @@ private: int *tag,*type,*mask,*image; double **x,**v,**f; + double **omega,**angmom; }; } diff -Naur lammps-11Nov09/src/comm.cpp lammps-12Nov09/src/comm.cpp --- lammps-11Nov09/src/comm.cpp 2009-09-22 17:08:40.000000000 -0600 +++ lammps-12Nov09/src/comm.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -63,6 +63,7 @@ multilo = multihi = NULL; cutghostmulti = NULL; cutghostuser = 0.0; + ghost_velocity = 0; // initialize comm buffers & exchange memory @@ -175,16 +176,28 @@ map_style = atom->map_style; // comm_only = 1 if only x,f are exchanged in forward/reverse comm + // comm_x_only = 0 if ghost_velocity since velocities are added comm_x_only = atom->avec->comm_x_only; comm_f_only = atom->avec->comm_f_only; + if (ghost_velocity) comm_x_only = 0; + + // set per-atom sizes for forward/reverse/border comm + // augment by velocity quantities if needed + + size_forward = atom->avec->size_forward; + size_reverse = atom->avec->size_reverse; + size_border = atom->avec->size_border; + + if (ghost_velocity) size_forward += atom->avec->size_velocity; + if (ghost_velocity) size_border += atom->avec->size_velocity; // maxforward = # of datums in largest forward communication // maxreverse = # of datums in largest reverse communication // query pair,fix,compute for their requirements - maxforward = MAX(atom->avec->size_comm,atom->avec->size_border); - maxreverse = atom->avec->size_reverse; + maxforward = MAX(size_forward,size_border); + maxreverse = size_reverse; if (force->pair) maxforward = MAX(maxforward,force->pair->comm_forward); if (force->pair) maxreverse = MAX(maxreverse,force->pair->comm_reverse); @@ -391,7 +404,7 @@ /* ---------------------------------------------------------------------- communication of atom coords every timestep - other stuff may also be sent via pack/unpack routines + other per-atom attributes may also be sent via pack/unpack routines ------------------------------------------------------------------------- */ void Comm::communicate() @@ -410,16 +423,24 @@ for (int iswap = 0; iswap < nswap; iswap++) { if (sendproc[iswap] != me) { if (comm_x_only) { - if (size_comm_recv[iswap]) buf = x[firstrecv[iswap]]; + if (size_forward_recv[iswap]) buf = x[firstrecv[iswap]]; else buf = NULL; - MPI_Irecv(buf,size_comm_recv[iswap],MPI_DOUBLE, + MPI_Irecv(buf,size_forward_recv[iswap],MPI_DOUBLE, recvproc[iswap],0,world,&request); n = avec->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); + } else if (ghost_velocity) { + MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE, + recvproc[iswap],0,world,&request); + n = avec->pack_comm_vel(sendnum[iswap],sendlist[iswap], + buf_send,pbc_flag[iswap],pbc[iswap]); + MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); + MPI_Wait(&request,&status); + avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_recv); } else { - MPI_Irecv(buf_recv,size_comm_recv[iswap],MPI_DOUBLE, + MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE, recvproc[iswap],0,world,&request); n = avec->pack_comm(sendnum[iswap],sendlist[iswap], buf_send,pbc_flag[iswap],pbc[iswap]); @@ -430,13 +451,19 @@ } else { if (comm_x_only) { - if (sendnum[iswap]) + if (sendnum[iswap]) { n = avec->pack_comm(sendnum[iswap],sendlist[iswap], - x[firstrecv[iswap]],pbc_flag[iswap],pbc[iswap]); - } else { - n = avec->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flag[iswap],pbc[iswap]); - avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send); + x[firstrecv[iswap]],pbc_flag[iswap], + pbc[iswap]); + } else if (ghost_velocity) { + n = avec->pack_comm_vel(sendnum[iswap],sendlist[iswap], + buf_send,pbc_flag[iswap],pbc[iswap]); + avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_send); + } else { + n = avec->pack_comm(sendnum[iswap],sendlist[iswap], + buf_send,pbc_flag[iswap],pbc[iswap]); + avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send); + } } } } @@ -444,7 +471,7 @@ /* ---------------------------------------------------------------------- reverse communication of forces on atoms every timestep - other stuff can also be sent via pack/unpack routines + other per-atom attributes may also be sent via pack/unpack routines ------------------------------------------------------------------------- */ void Comm::reverse_communicate() @@ -623,7 +650,6 @@ MPI_Request request; MPI_Status status; AtomVec *avec = atom->avec; - int size_border = avec->size_border; // clear old ghosts @@ -719,8 +745,12 @@ if (nsend*size_border > maxsend) grow_send(nsend*size_border,0); - n = avec->pack_border(nsend,sendlist[iswap],buf_send, - pbc_flag[iswap],pbc[iswap]); + if (ghost_velocity) + n = avec->pack_border_vel(nsend,sendlist[iswap],buf_send, + pbc_flag[iswap],pbc[iswap]); + else + n = avec->pack_border(nsend,sendlist[iswap],buf_send, + pbc_flag[iswap],pbc[iswap]); // swap atoms with other proc // put incoming ghosts at end of my atom arrays @@ -743,7 +773,10 @@ // unpack buffer - avec->unpack_border(nrecv,atom->nlocal+atom->nghost,buf); + if (ghost_velocity) + avec->unpack_border_vel(nrecv,atom->nlocal+atom->nghost,buf); + else + avec->unpack_border(nrecv,atom->nlocal+atom->nghost,buf); // set all pointers & counters @@ -751,9 +784,9 @@ rmax = MAX(rmax,nrecv); sendnum[iswap] = nsend; recvnum[iswap] = nrecv; - size_comm_recv[iswap] = nrecv * avec->size_comm; - size_reverse_send[iswap] = nrecv * avec->size_reverse; - size_reverse_recv[iswap] = nsend * avec->size_reverse; + size_forward_recv[iswap] = nrecv*size_forward; + size_reverse_send[iswap] = nrecv*size_reverse; + size_reverse_recv[iswap] = nsend*size_reverse; firstrecv[iswap] = atom->nlocal + atom->nghost; atom->nghost += nrecv; iswap++; @@ -1469,7 +1502,7 @@ recvnum = (int *) memory->smalloc(n*sizeof(int),"comm:recvnum"); sendproc = (int *) memory->smalloc(n*sizeof(int),"comm:sendproc"); recvproc = (int *) memory->smalloc(n*sizeof(int),"comm:recvproc"); - size_comm_recv = (int *) memory->smalloc(n*sizeof(int),"comm:size"); + size_forward_recv = (int *) memory->smalloc(n*sizeof(int),"comm:size"); size_reverse_send = (int *) memory->smalloc(n*sizeof(int),"comm:size"); size_reverse_recv = (int *) memory->smalloc(n*sizeof(int),"comm:size"); slablo = (double *) memory->smalloc(n*sizeof(double),"comm:slablo"); @@ -1500,7 +1533,7 @@ memory->sfree(recvnum); memory->sfree(sendproc); memory->sfree(recvproc); - memory->sfree(size_comm_recv); + memory->sfree(size_forward_recv); memory->sfree(size_reverse_send); memory->sfree(size_reverse_recv); memory->sfree(slablo); @@ -1549,6 +1582,12 @@ if (cutghostuser < 0.0) error->all("Invalid cutoff in communicate command"); iarg += 2; + } else if (strcmp(arg[iarg],"vel") == 0) { + if (iarg+2 > narg) error->all("Illegal communicate command"); + if (strcmp(arg[iarg+1],"yes") == 0) ghost_velocity = 1; + else if (strcmp(arg[iarg+1],"yes") == 0) ghost_velocity = 0; + else error->all("Illegal communicate command"); + iarg += 2; } else error->all("Illegal communicate command"); } } diff -Naur lammps-11Nov09/src/comm.h lammps-12Nov09/src/comm.h --- lammps-11Nov09/src/comm.h 2008-12-02 13:18:37.000000000 -0700 +++ lammps-12Nov09/src/comm.h 2009-11-09 11:20:20.000000000 -0700 @@ -28,12 +28,15 @@ int procneigh[3][2]; // my 6 neighboring procs int nswap; // # of swaps to perform int need[3]; // procs I need atoms from in each dim + int ghost_velocity; // 1 if ghost atoms have velocity, 0 if not int maxforward_fix; // comm sizes called from Fix,Pair int maxreverse_fix; int maxforward_pair; int maxreverse_pair; double cutghost[3]; // cutoffs used for acquiring ghost atoms + double cutghostuser; // user-specified ghost cutoff + Comm(class LAMMPS *); ~Comm(); @@ -60,9 +63,12 @@ private: int triclinic; // 0 if domain is orthog, 1 if triclinic int maxswap; // max # of swaps memory is allocated for + int size_forward; // # of per-atom datums in forward comm + int size_reverse; // # of datums in reverse comm + int size_border; // # of datums in forward border comm int *sendnum,*recvnum; // # of atoms to send/recv in each swap int *sendproc,*recvproc; // proc to send/recv to/from at each swap - int *size_comm_recv; // # of values to recv in each forward comm + int *size_forward_recv; // # of values to recv in each forward comm int *size_reverse_send; // # to send in each reverse comm int *size_reverse_recv; // # to recv in each reverse comm double *slablo,*slabhi; // bounds of slab to send at each swap @@ -74,7 +80,6 @@ int map_style; // non-0 if global->local mapping is done int ***grid2proc; // which proc owns i,j,k loc in 3d grid int bordergroup; // only communicate this group in borders - double cutghostuser; // user-specified ghost cutoff int *firstrecv; // where to put 1st recv atom in each swap int **sendlist; // list of atoms to send in each swap diff -Naur lammps-11Nov09/src/compute_heat_flux.cpp lammps-12Nov09/src/compute_heat_flux.cpp --- lammps-11Nov09/src/compute_heat_flux.cpp 2009-07-02 12:01:55.000000000 -0600 +++ lammps-12Nov09/src/compute_heat_flux.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -24,6 +24,7 @@ #include "atom_vec.h" #include "update.h" #include "force.h" +#include "comm.h" #include "pair.h" #include "modify.h" #include "group.h" @@ -75,7 +76,7 @@ { // error checks - if (atom->avec->ghost_velocity == 0) + if (comm->ghost_velocity == 0) error->all("Compute heat/flux requires ghost atoms store velocity"); if (force->pair == NULL || force->pair->single_enable == 0) diff -Naur lammps-11Nov09/src/fix_langevin.cpp lammps-12Nov09/src/fix_langevin.cpp --- lammps-11Nov09/src/fix_langevin.cpp 2009-08-18 12:11:26.000000000 -0600 +++ lammps-12Nov09/src/fix_langevin.cpp 2009-11-06 16:34:26.000000000 -0700 @@ -450,11 +450,11 @@ strcpy(id_temp,arg[1]); int icompute = modify->find_compute(id_temp); - if (icompute < 0) error->all("Could not find fix_modify temp ID"); + if (icompute < 0) error->all("Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; if (temperature->tempflag == 0) - error->all("Fix_modify temp ID does not compute temperature"); + error->all("Fix_modify temperature ID does not compute temperature"); if (temperature->igroup != igroup && comm->me == 0) error->warning("Group for fix_modify temp != fix group"); return 2; diff -Naur lammps-11Nov09/src/pair_dpd.cpp lammps-12Nov09/src/pair_dpd.cpp --- lammps-11Nov09/src/pair_dpd.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-12Nov09/src/pair_dpd.cpp 2009-11-09 11:20:20.000000000 -0700 @@ -0,0 +1,388 @@ +/* ---------------------------------------------------------------------- + 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: Kurt Smith (U Pittsburgh) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "pair_dpd.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "update.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "random_mars.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +#define EPSILON 1.0e-10 + +/* ---------------------------------------------------------------------- */ + +PairDPD::PairDPD(LAMMPS *lmp) : Pair(lmp) +{ + random = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairDPD::~PairDPD() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + + memory->destroy_2d_double_array(cut); + memory->destroy_2d_double_array(a0); + memory->destroy_2d_double_array(gamma); + memory->destroy_2d_double_array(sigma); + } + + if (random) delete random; +} + +/* ---------------------------------------------------------------------- */ + +void PairDPD::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double vxtmp,vytmp,vztmp,delvx,delvy,delvz; + double rsq,r,rinv,dot,wd,randnum,factor_dpd; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double dtinvsqrt = 1.0/sqrt(update->dt); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + if (j < nall) factor_dpd = 1.0; + else { + factor_dpd = special_lj[j/nall]; + j %= nall; + } + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + if (r < EPSILON) continue; // r can be 0.0 in DPD systems + rinv = 1.0/r; + delvx = vxtmp - v[j][0]; + delvy = vytmp - v[j][1]; + delvz = vztmp - v[j][2]; + dot = delx*delvx + dely*delvy + delz*delvz; + wd = 1.0 - r/cut[itype][jtype]; + randnum = random->gaussian(); + + // conservative force = a0 * wd + // drag force = -gamma * wd^2 * (delx dot delv) / r + // random force = sigma * wd * rnd * dtinvsqrt; + + fpair = a0[itype][jtype]*wd; + fpair -= gamma[itype][jtype]*wd*wd*dot*rinv; + fpair += sigma[itype][jtype]*wd*randnum*dtinvsqrt; + fpair *= factor_dpd*rinv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + evdwl = -a0[itype][jtype] * r * (1.0 - 0.5*r/cut[itype][jtype]); + evdwl *= factor_dpd; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairDPD::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); + a0 = memory->create_2d_double_array(n+1,n+1,"pair:a0"); + gamma = memory->create_2d_double_array(n+1,n+1,"pair:gamma"); + sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairDPD::settings(int narg, char **arg) +{ + if (narg != 3) error->all("Illegal pair_style command"); + + temperature = force->numeric(arg[0]); + cut_global = force->numeric(arg[1]); + seed = force->inumeric(arg[2]); + + // initialize Marsaglia RNG with processor-unique seed + + if (seed <= 0) error->all("Illegal pair_style command"); + delete random; + random = new RanMars(lmp,seed + comm->me); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairDPD::coeff(int narg, char **arg) +{ + if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double a0_one = force->numeric(arg[2]); + double gamma_one = force->numeric(arg[3]); + + double cut_one = cut_global; + if (narg == 5) cut_one = force->numeric(arg[4]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + a0[i][j] = a0_one; + gamma[i][j] = gamma_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairDPD::init_style() +{ + if (comm->ghost_velocity == 0) + error->all("Pair dpd requires ghost atoms store velocity"); + + // if newton off, forces between atoms ij will be double computed + // using different random numbers + + if (force->newton_pair == 0 && comm->me == 0) error->warning( + "Pair dpd needs newton pair on for momentum conservation"); + + int irequest = neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairDPD::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); + + sigma[i][j] = sqrt(2.0*force->boltz*temperature*gamma[i][j]); + + cut[j][i] = cut[i][j]; + a0[j][i] = a0[i][j]; + gamma[j][i] = gamma[i][j]; + sigma[j][i] = sigma[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDPD::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&a0[i][j],sizeof(double),1,fp); + fwrite(&gamma[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDPD::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&a0[i][j],sizeof(double),1,fp); + fread(&gamma[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&gamma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDPD::write_restart_settings(FILE *fp) +{ + fwrite(&temperature,sizeof(double),1,fp); + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&seed,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDPD::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&temperature,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&seed,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&temperature,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&seed,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + + // initialize Marsaglia RNG with processor-unique seed + // same seed that pair_style command initially specified + + if (random) delete random; + random = new RanMars(lmp,seed + comm->me); +} + +/* ---------------------------------------------------------------------- */ + +double PairDPD::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_dpd, double &fforce) +{ + double r,rinv,wd,phi; + + r = sqrt(rsq); + if (r < EPSILON) { + fforce = 0.0; + return 0.0; + } + + rinv = 1.0/r; + wd = 1.0 - r/cut[itype][jtype]; + fforce = a0[itype][jtype]*wd * factor_dpd*rinv; + + phi = -a0[itype][jtype] * r * (1.0 - 0.5*r/cut[itype][jtype]); + return factor_dpd*phi; +} diff -Naur lammps-11Nov09/src/pair_dpd.h lammps-12Nov09/src/pair_dpd.h --- lammps-11Nov09/src/pair_dpd.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-12Nov09/src/pair_dpd.h 2009-11-09 11:20:20.000000000 -0700 @@ -0,0 +1,49 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_DPD_H +#define PAIR_DPD_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairDPD : public Pair { + public: + PairDPD(class LAMMPS *); + ~PairDPD(); + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + double single(int, int, int, int, double, double, double, double &); + + private: + double cut_global,temperature; + int seed; + double **cut; + double **a0,**gamma; + double **sigma; + class RanMars *random; + + void allocate(); +}; + +} + +#endif diff -Naur lammps-11Nov09/src/style.h lammps-12Nov09/src/style.h --- lammps-11Nov09/src/style.h 2009-11-06 09:48:55.000000000 -0700 +++ lammps-12Nov09/src/style.h 2009-11-09 11:20:20.000000000 -0700 @@ -312,6 +312,7 @@ #include "pair_buck_coul_cut.h" #include "pair_coul_cut.h" #include "pair_coul_debye.h" +#include "pair_dpd.h" #include "pair_hybrid.h" #include "pair_hybrid_overlay.h" #include "pair_lj_cut.h" @@ -333,6 +334,7 @@ PairStyle(buck/coul/cut,PairBuckCoulCut) PairStyle(coul/cut,PairCoulCut) PairStyle(coul/debye,PairCoulDebye) +PairStyle(dpd,PairDPD) PairStyle(hybrid,PairHybrid) PairStyle(hybrid/overlay,PairHybridOverlay) PairStyle(lj/cut,PairLJCut) diff -Naur lammps-11Nov09/src/style_dpd.h lammps-12Nov09/src/style_dpd.h --- lammps-11Nov09/src/style_dpd.h 2009-11-09 11:42:10.000000000 -0700 +++ lammps-12Nov09/src/style_dpd.h 2009-09-28 10:08:11.000000000 -0600 @@ -0,0 +1,28 @@ +/* ---------------------------------------------------------------------- + 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 AtomInclude +#include "atom_vec_dpd.h" +#endif + +#ifdef AtomClass +AtomStyle(dpd,AtomVecDPD) +#endif + +#ifdef PairInclude +#include "pair_dpd.h" +#endif + +#ifdef PairClass +PairStyle(dpd,PairDPD) +#endif diff -Naur lammps-11Nov09/src/style_packages.h lammps-12Nov09/src/style_packages.h --- lammps-11Nov09/src/style_packages.h 2009-11-06 09:48:55.000000000 -0700 +++ lammps-12Nov09/src/style_packages.h 2009-11-09 11:20:20.000000000 -0700 @@ -17,7 +17,6 @@ #include "style_class2.h" #include "style_colloid.h" #include "style_dipole.h" -#include "style_dpd.h" #include "style_dsmc.h" #include "style_gpu.h" #include "style_granular.h"