added false positive, removed graphs from examples, updated langevin kokkos, improved diff readability in langevin

This commit is contained in:
charlie sievers 2019-09-12 16:34:15 -07:00
parent 1e8e34f33d
commit a948a34f8a
15 changed files with 420 additions and 191 deletions

3
.gitignore vendored
View File

@ -43,6 +43,3 @@ Thumbs.db
/Makefile
/cmake_install.cmake
/lmp
#python example
/example/python/gjf_python

View File

@ -2269,6 +2269,7 @@ qoffload
qopenmp
qoverride
qtb
quadratically
quadrupolar
Quant
quartic

View File

@ -1,5 +0,0 @@
# LAMMPS GJF-2GJ MOLECULAR DYNAMICS RESULTS
## GJF-2GJ THERMOSTAT
This directory contains a series of graphs, which display the results from numermous molecular dynamics simulations. All simulations are run in the NVT ensemble. Two systems are reported, guaiacol and argon. The damping parameter and the timestep are varied. Also the temperature is varied for argon. GJF U is the half-step velocity (vhalf) and GJF V is the onsite velocity (vfull). GJF U and GJF V represent exactly the same configurational statistics.

View File

@ -11,20 +11,23 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_langevin_kokkos.h"
#include <cmath>
#include <cstdio>
#include <cstring>
#include "fix_langevin_kokkos.h"
#include "atom_masks.h"
#include "atom_kokkos.h"
#include "force.h"
#include "group.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "memory_kokkos.h"
#include "group.h"
#include "random_mars.h"
#include "compute.h"
#include "comm.h"
#include "modify.h"
#include "input.h"
#include "region.h"
#include "variable.h"
using namespace LAMMPS_NS;
@ -61,7 +64,6 @@ FixLangevinKokkos<DeviceType>::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a
k_ratio.template modify<LMPHostType>();
if(gjfflag){
nvalues = 3;
grow_arrays(atomKK->nmax);
atom->add_callback(0);
// initialize franprev to zero
@ -69,8 +71,12 @@ FixLangevinKokkos<DeviceType>::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a
franprev[i][0] = 0.0;
franprev[i][1] = 0.0;
franprev[i][2] = 0.0;
lv[i][0] = 0.0;
lv[i][1] = 0.0;
lv[i][2] = 0.0;
}
k_franprev.template modify<LMPHostType>();
k_lv.template modify<LMPHostType>();
}
if(zeroflag){
k_fsumall = tdual_double_1d_3n("langevin:fsumall");
@ -94,6 +100,7 @@ FixLangevinKokkos<DeviceType>::~FixLangevinKokkos()
memoryKK->destroy_kokkos(k_ratio,ratio);
memoryKK->destroy_kokkos(k_flangevin,flangevin);
if(gjfflag) memoryKK->destroy_kokkos(k_franprev,franprev);
if(gjfflag) memoryKK->destroy_kokkos(k_lv,lv);
memoryKK->destroy_kokkos(k_tforce,tforce);
}
@ -107,6 +114,11 @@ void FixLangevinKokkos<DeviceType>::init()
error->all(FLERR,"Fix langevin omega is not yet implemented with kokkos");
if(ascale)
error->all(FLERR,"Fix langevin angmom is not yet implemented with kokkos");
if(gjfflag && tbiasflag)
error->all(FLERR,"Fix langevin gjf + tbias is not yet implemented with kokkos");
if(gjfflag && tbiasflag)
error->warning(FLERR,"Fix langevin gjf + kokkos is not implemented with random gaussians,"
" this may cause errors in kinetic fluctuations");
// prefactors are modified in the init
k_gfactor1.template modify<LMPHostType>();
@ -121,6 +133,42 @@ void FixLangevinKokkos<DeviceType>::grow_arrays(int nmax)
memoryKK->grow_kokkos(k_franprev,franprev,nmax,3,"langevin:franprev");
d_franprev = k_franprev.template view<DeviceType>();
h_franprev = k_franprev.template view<LMPHostType>();
memoryKK->grow_kokkos(k_lv,lv,nmax,3,"langevin:lv");
d_lv = k_lv.template view<DeviceType>();
h_lv = k_lv.template view<LMPHostType>();
}
/* ----------------------------------------------------------------------
allow for both per-type and per-atom mass
------------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::initial_integrate(int vflag)
{
atomKK->sync(execution_space,datamask_read);
atomKK->modified(execution_space,datamask_modify);
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
int nlocal = atomKK->nlocal;
if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
FixLangevinKokkosInitialIntegrateFunctor<DeviceType> functor(this);
Kokkos::parallel_for(nlocal,functor);
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::initial_integrate_item(int i) const
{
if (mask[i] & groupbit) {
f(i,0) /= gjfa;
f(i,1) /= gjfa;
f(i,2) /= gjfa;
v(i,0) = d_lv(i,0);
v(i,1) = d_lv(i,1);
v(i,2) = d_lv(i,2);
}
}
/* ---------------------------------------------------------------------- */
@ -140,6 +188,7 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
k_gfactor2.template sync<DeviceType>();
k_ratio.template sync<DeviceType>();
if(gjfflag) k_franprev.template sync<DeviceType>();
if(gjfflag) k_lv.template sync<DeviceType>();
boltz = force->boltz;
dt = update->dt;
@ -177,7 +226,7 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
atomKK->sync(temperature->execution_space,temperature->datamask_read);
temperature->compute_scalar();
temperature->remove_bias_all(); // modifies velocities
// if temeprature compute is kokkosized host-devcie comm won't be needed
// if temeprature compute is kokkosized host-device comm won't be needed
atomKK->modified(temperature->execution_space,temperature->datamask_modify);
atomKK->sync(execution_space,temperature->datamask_modify);
}
@ -481,6 +530,7 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
// set modify flags for the views modified in post_force functor
if (gjfflag) k_franprev.template modify<DeviceType>();
if (gjfflag) k_lv.template modify<DeviceType>();
if (tallyflag) k_flangevin.template modify<DeviceType>();
// set total force to zero
@ -550,6 +600,10 @@ FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
}
if (Tp_GJF) {
d_lv(i,0) = gjfsib*v(i,0);
d_lv(i,1) = gjfsib*v(i,1);
d_lv(i,2) = gjfsib*v(i,2);
fswap = 0.5*(fran[0]+d_franprev(i,0));
d_franprev(i,0) = fran[0];
fran[0] = fswap;
@ -560,15 +614,15 @@ FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
d_franprev(i,2) = fran[2];
fran[2] = fswap;
fdrag[0] *= gjffac;
fdrag[1] *= gjffac;
fdrag[2] *= gjffac;
fran[0] *= gjffac;
fran[1] *= gjffac;
fran[2] *= gjffac;
f(i,0) *= gjffac;
f(i,1) *= gjffac;
f(i,2) *= gjffac;
fdrag[0] *= gjfa;
fdrag[1] *= gjfa;
fdrag[2] *= gjfa;
fran[0] *= gjfa;
fran[1] *= gjfa;
fran[2] *= gjfa;
f(i,0) *= gjfa;
f(i,1) *= gjfa;
f(i,2) *= gjfa;
}
f(i,0) += fdrag[0] + fran[0];
@ -576,6 +630,17 @@ FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
f(i,2) += fdrag[2] + fran[2];
if (Tp_TALLY) {
if (Tp_GJF){
fdrag[0] = gamma1*d_lv(i,0)/gjfsib/gjfsib;
fdrag[1] = gamma1*d_lv(i,1)/gjfsib/gjfsib;
fdrag[2] = gamma1*d_lv(i,2)/gjfsib/gjfsib;
fswap = (2*fran[0]/gjfa - d_franprev(i,0))/gjfsib;
fran[0] = fswap;
fswap = (2*fran[1]/gjfa - d_franprev(i,1))/gjfsib;
fran[1] = fswap;
fswap = (2*fran[2]/gjfa - d_franprev(i,2))/gjfsib;
fran[2] = fswap;
}
d_flangevin(i,0) = fdrag[0] + fran[0];
d_flangevin(i,1) = fdrag[1] + fran[1];
d_flangevin(i,2) = fdrag[2] + fran[2];
@ -719,9 +784,10 @@ double FixLangevinKokkos<DeviceType>::compute_energy_item(int i) const
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::end_of_step()
{
if (!tallyflag) return;
if (!tallyflag && !gjfflag) return;
v = atomKK->k_v.template view<DeviceType>();
f = atomKK->k_f.template view<DeviceType>();
mask = atomKK->k_mask.template view<DeviceType>();
atomKK->sync(execution_space,V_MASK | MASK_MASK);
@ -733,9 +799,81 @@ void FixLangevinKokkos<DeviceType>::end_of_step()
FixLangevinKokkosTallyEnergyFunctor<DeviceType> tally_functor(this);
Kokkos::parallel_reduce(nlocal,tally_functor,energy_onestep);
if (gjfflag){
if (rmass.data()) {
FixLangevinKokkosEndOfStepFunctor<DeviceType,1> functor(this);
Kokkos::parallel_for(nlocal,functor);
} else {
mass = atomKK->k_mass.view<DeviceType>();
FixLangevinKokkosEndOfStepFunctor<DeviceType,0> functor(this);
Kokkos::parallel_for(nlocal,functor);
}
}
energy += energy_onestep*update->dt;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::end_of_step_item(int i) const {
double tmp[3];
if (mask[i] & groupbit) {
const double dtfm = force->ftm2v * 0.5 * dt / mass[type[i]];
tmp[0] = v(i,0);
tmp[1] = v(i,1);
tmp[2] = v(i,2);
if (!fsflag){
v(i,0) = d_lv(i,0);
v(i,1) = d_lv(i,1);
v(i,2) = d_lv(i,2);
} else {
v(i,0) = 0.5 * gjfsib * gjfsib * (v(i,0) + dtfm * f(i,0) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,0)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,0);
v(i,1) = 0.5 * gjfsib * gjfsib * (v(i,1) + dtfm * f(i,1) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,1)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,1);
v(i,2) = 0.5 * gjfsib * gjfsib * (v(i,2) + dtfm * f(i,2) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,2)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,2);
}
d_lv(i,0) = tmp[0];
d_lv(i,1) = tmp[1];
d_lv(i,2) = tmp[2];
}
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::end_of_step_rmass_item(int i) const
{
double tmp[3];
if (mask[i] & groupbit) {
const double dtfm = force->ftm2v * 0.5 * dt / rmass[i];
tmp[0] = v(i,0);
tmp[1] = v(i,1);
tmp[2] = v(i,2);
if (!fsflag){
v(i,0) = d_lv(i,0);
v(i,1) = d_lv(i,1);
v(i,2) = d_lv(i,2);
} else {
v(i,0) = 0.5 * gjfsib * gjfsib * (v(i,0) + dtfm * f(i,0) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,0)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,0);
v(i,1) = 0.5 * gjfsib * gjfsib * (v(i,1) + dtfm * f(i,1) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,1) - d_franprev(i,1)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,1);
v(i,2) = 0.5 * gjfsib * gjfsib * (v(i,2) + dtfm * f(i,2) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,2) - d_franprev(i,2)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,2);
}
d_lv(i,0) = tmp[0];
d_lv(i,1) = tmp[1];
d_lv(i,2) = tmp[2];
}
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
@ -743,10 +881,15 @@ void FixLangevinKokkos<DeviceType>::end_of_step()
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < nvalues; m++)
h_franprev(j,m) = h_franprev(i,m);
h_franprev(j,0) = h_franprev(i,0);
h_franprev(j,1) = h_franprev(i,1);
h_franprev(j,2) = h_franprev(i,2);
h_lv(j,0) = h_lv(i,0);
h_lv(j,1) = h_lv(i,1);
h_lv(j,2) = h_lv(i,2);
k_franprev.template modify<LMPHostType>();
k_lv.template modify<LMPHostType>();
}
@ -765,6 +908,7 @@ void FixLangevinKokkos<DeviceType>::cleanup_copy()
tforce = NULL;
gjfflag = 0;
franprev = NULL;
lv = NULL;
id = style = NULL;
vatom = NULL;
}

View File

@ -56,6 +56,9 @@ namespace LAMMPS_NS {
template<class DeviceType>
class FixLangevinKokkos;
template <class DeviceType>
class FixLangevinKokkosInitialIntegrateFunctor;
template<class DeviceType,int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
class FixLangevinKokkosPostForceFunctor;
@ -72,6 +75,7 @@ namespace LAMMPS_NS {
void cleanup_copy();
void init();
void initial_integrate(int);
void post_force(int);
void reset_dt();
void grow_arrays(int);
@ -79,6 +83,12 @@ namespace LAMMPS_NS {
double compute_scalar();
void end_of_step();
KOKKOS_INLINE_FUNCTION
void initial_integrate_item(int) const;
KOKKOS_INLINE_FUNCTION
void initial_integrate_rmass_item(int) const;
template<int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
KOKKOS_INLINE_FUNCTION
@ -90,14 +100,25 @@ namespace LAMMPS_NS {
KOKKOS_INLINE_FUNCTION
double compute_energy_item(int) const;
KOKKOS_INLINE_FUNCTION
void end_of_step_item(int) const;
KOKKOS_INLINE_FUNCTION
void end_of_step_rmass_item(int) const;
private:
class CommKokkos *commKK;
typename ArrayTypes<DeviceType>::t_float_1d rmass;
typename ArrayTypes<DeviceType>::t_float_1d mass;
typename ArrayTypes<DeviceType>::tdual_double_2d k_franprev;
typename ArrayTypes<DeviceType>::t_double_2d d_franprev;
HAT::t_double_2d h_franprev;
typename ArrayTypes<DeviceType>::tdual_double_2d k_lv;
typename ArrayTypes<DeviceType>::t_double_2d d_lv;
HAT::t_double_2d h_lv;
typename ArrayTypes<DeviceType>::tdual_double_2d k_flangevin;
typename ArrayTypes<DeviceType>::t_double_2d d_flangevin;
HAT::t_double_2d h_flangevin;
@ -130,6 +151,21 @@ namespace LAMMPS_NS {
};
template <class DeviceType>
struct FixLangevinKokkosInitialIntegrateFunctor {
typedef DeviceType device_type ;
FixLangevinKokkos<DeviceType> c;
FixLangevinKokkosInitialIntegrateFunctor(FixLangevinKokkos<DeviceType>* c_ptr):
c(*c_ptr) {c.cleanup_copy();};
KOKKOS_INLINE_FUNCTION
void operator()(const int i) const {
c.initial_integrate_item(i);
}
};
template <class DeviceType,int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
struct FixLangevinKokkosPostForceFunctor {
@ -207,6 +243,21 @@ namespace LAMMPS_NS {
update += source;
}
};
template <class DeviceType, int RMass>
struct FixLangevinKokkosEndOfStepFunctor {
typedef DeviceType device_type ;
FixLangevinKokkos<DeviceType> c;
FixLangevinKokkosEndOfStepFunctor(FixLangevinKokkos<DeviceType>* c_ptr):
c(*c_ptr) {c.cleanup_copy();}
KOKKOS_INLINE_FUNCTION
void operator()(const int i) const {
if (RMass) c.end_of_step_rmass_item(i);
else c.end_of_step_item(i);
}
};
}
#endif

View File

@ -223,7 +223,7 @@ void FixLangevin::init()
{
if (gjfflag){
if (t_period*2 == update->dt)
error->all(FLERR,"Fix langevin gjf cannot have t_period equal to dt/2 at the start");
error->all(FLERR,"Fix langevin gjf cannot have t_period equal to dt/2");
// warn if any integrate fix comes after this one
int before = 1;
@ -440,124 +440,124 @@ void FixLangevin::post_force(int /*vflag*/)
if (zeroflag) post_force_templated<1,1,1,1,1,1>();
else post_force_templated<1,1,1,1,1,0>();
else
if (zeroflag) post_force_templated<1,1,1,1,0,1>();
else post_force_templated<1,1,1,1,0,0>();
if (zeroflag) post_force_templated<1,1,1,1,0,1>();
else post_force_templated<1,1,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,1,1,0,1,1>();
else post_force_templated<1,1,1,0,1,0>();
if (rmass)
if (zeroflag) post_force_templated<1,1,1,0,1,1>();
else post_force_templated<1,1,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,1,0,0,1>();
else post_force_templated<1,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,0,1,1,1>();
else post_force_templated<1,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,1,0,1>();
else post_force_templated<1,1,0,1,0,0>();
else
if (zeroflag) post_force_templated<1,1,1,0,0,1>();
else post_force_templated<1,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,0,1,1,1>();
else post_force_templated<1,1,0,1,1,0>();
if (rmass)
if (zeroflag) post_force_templated<1,1,0,0,1,1>();
else post_force_templated<1,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,0,0,1>();
else post_force_templated<1,1,0,0,0,0>();
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,1,1,1,1>();
else post_force_templated<1,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,1,0,1>();
else post_force_templated<1,0,1,1,0,0>();
else
if (zeroflag) post_force_templated<1,1,0,1,0,1>();
else post_force_templated<1,1,0,1,0,0>();
if (rmass)
if (zeroflag) post_force_templated<1,0,1,0,1,1>();
else post_force_templated<1,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,0,0,1>();
else post_force_templated<1,0,1,0,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,1,0,0,1,1>();
else post_force_templated<1,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,0,0,1>();
else post_force_templated<1,1,0,0,0,0>();
else
if (tallyflag || fsflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,1,1,1,1>();
else post_force_templated<1,0,1,1,1,0>();
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,0,1,1,1>();
else post_force_templated<1,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,1,0,1>();
else post_force_templated<1,0,0,1,0,0>();
else
if (zeroflag) post_force_templated<1,0,1,1,0,1>();
else post_force_templated<1,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,0,1,0,1,1>();
else post_force_templated<1,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,0,0,1>();
else post_force_templated<1,0,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,0,1,1,1>();
else post_force_templated<1,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,1,0,1>();
else post_force_templated<1,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,0,0,0,1,1>();
else post_force_templated<1,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,0,0,1>();
else post_force_templated<1,0,0,0,0,0>();
if (rmass)
if (zeroflag) post_force_templated<1,0,0,0,1,1>();
else post_force_templated<1,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,0,0,1>();
else post_force_templated<1,0,0,0,0,0>();
else
if (gjfflag)
if (tallyflag || fsflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,1,1,1,1>();
else post_force_templated<0,1,1,1,1,0>();
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,1,1,1,1>();
else post_force_templated<0,1,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,1,1,0,1>();
else post_force_templated<0,1,1,1,0,0>();
else
if (zeroflag) post_force_templated<0,1,1,1,0,1>();
else post_force_templated<0,1,1,1,0,0>();
if (rmass)
if (zeroflag) post_force_templated<0,1,1,0,1,1>();
else post_force_templated<0,1,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,1,0,0,1>();
else post_force_templated<0,1,1,0,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,1,0,1,1>();
else post_force_templated<0,1,1,0,1,0>();
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,0,1,1,1>();
else post_force_templated<0,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,1,0,1>();
else post_force_templated<0,1,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,0,0,1,1>();
else post_force_templated<0,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,0,0,1>();
else post_force_templated<0,1,0,0,0,0>();
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,1,1,1,1>();
else post_force_templated<0,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,1,0,1>();
else post_force_templated<0,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,1,0,1,1>();
else post_force_templated<0,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,0,0,1>();
else post_force_templated<0,0,1,0,0,0>();
else
if (zeroflag) post_force_templated<0,1,1,0,0,1>();
else post_force_templated<0,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,0,1,1,1>();
else post_force_templated<0,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,1,0,1>();
else post_force_templated<0,1,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,0,0,1,1>();
else post_force_templated<0,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,0,0,1>();
else post_force_templated<0,1,0,0,0,0>();
else
if (tallyflag || fsflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,1,1,1,1>();
else post_force_templated<0,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,1,0,1>();
else post_force_templated<0,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,1,0,1,1>();
else post_force_templated<0,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,0,0,1>();
else post_force_templated<0,0,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,0,1,1,1>();
else post_force_templated<0,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,1,0,1>();
else post_force_templated<0,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,0,0,1,1>();
else post_force_templated<0,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,0,0,1>();
else post_force_templated<0,0,0,0,0,0>();
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,0,1,1,1>();
else post_force_templated<0,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,1,0,1>();
else post_force_templated<0,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,0,0,1,1>();
else post_force_templated<0,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,0,0,1>();
else post_force_templated<0,0,0,0,0,0>();
}
/* ---------------------------------------------------------------------- */
@ -572,7 +572,7 @@ void FixLangevin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
------------------------------------------------------------------------- */
template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
void FixLangevin::post_force_templated()
{
double gamma1,gamma2;
@ -802,9 +802,9 @@ void FixLangevin::compute_target()
input->variable->compute_atom(tvar,igroup,tforce,1,0);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (tforce[i] < 0.0)
error->one(FLERR,
"Fix langevin variable returned negative temperature");
if (tforce[i] < 0.0)
error->one(FLERR,
"Fix langevin variable returned negative temperature");
}
modify->addstep_compute(update->ntimestep + 1);
}
@ -1001,8 +1001,8 @@ void FixLangevin::reset_dt()
if (atom->mass) {
for (int i = 1; i <= atom->ntypes; i++) {
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
gfactor2[i] *= 1.0/sqrt(ratio[i]);
}
}

View File

@ -2,10 +2,12 @@
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.
------------------------------------------------------------------------- */
@ -22,63 +24,63 @@ FixStyle(langevin,FixLangevin)
namespace LAMMPS_NS {
class FixLangevin : public Fix {
public:
FixLangevin(class LAMMPS *, int, char **);
virtual ~FixLangevin();
int setmask();
void init();
void setup(int);
virtual void initial_integrate(int);
virtual void post_force(int);
void post_force_respa(int, int, int);
virtual void end_of_step();
void reset_target(double);
void reset_dt();
int modify_param(int, char **);
virtual double compute_scalar();
double memory_usage();
virtual void *extract(const char *, int &);
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
class FixLangevin : public Fix {
public:
FixLangevin(class LAMMPS *, int, char **);
virtual ~FixLangevin();
int setmask();
void init();
void setup(int);
virtual void initial_integrate(int);
virtual void post_force(int);
void post_force_respa(int, int, int);
virtual void end_of_step();
void reset_target(double);
void reset_dt();
int modify_param(int, char **);
virtual double compute_scalar();
double memory_usage();
virtual void *extract(const char *, int &);
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
protected:
int gjfflag,fsflag,oflag,tallyflag,zeroflag,tbiasflag;
int flangevin_allocated;
double ascale;
double t_start,t_stop,t_period,t_target;
double *gfactor1,*gfactor2,*ratio;
double energy,energy_onestep;
double tsqrt;
int tstyle,tvar;
double gjfa, gjfsib; //gjf a and gjf sqrt inverse b
char *tstr;
protected:
int gjfflag,fsflag,oflag,tallyflag,zeroflag,tbiasflag;
int flangevin_allocated;
double ascale;
double t_start,t_stop,t_period,t_target;
double *gfactor1,*gfactor2,*ratio;
double energy,energy_onestep;
double tsqrt;
int tstyle,tvar;
double gjfa, gjfsib; //gjf a and gjf sqrt inverse b
char *tstr;
class AtomVecEllipsoid *avec;
class AtomVecEllipsoid *avec;
int maxatom1,maxatom2;
double **flangevin;
double *tforce;
double **franprev;
double **lv; //half step velocity
int maxatom1,maxatom2;
double **flangevin;
double *tforce;
double **franprev;
double **lv; //half step velocity
char *id_temp;
class Compute *temperature;
char *id_temp;
class Compute *temperature;
int nlevels_respa;
class RanMars *random;
int seed;
int nlevels_respa;
class RanMars *random;
int seed;
template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
void post_force_templated();
template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
void post_force_templated();
void omega_thermostat();
void angmom_thermostat();
void compute_target();
};
void omega_thermostat();
void angmom_thermostat();
void compute_target();
};
}
@ -86,35 +88,74 @@ namespace LAMMPS_NS {
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix langevin period must be > 0.0
The time window for temperature relaxation must be > 0
E: Fix langevin omega requires atom style sphere
Self-explanatory.
E: Fix langevin angmom requires atom style ellipsoid
Self-explanatory.
E: Variable name for fix langevin does not exist
Self-explanatory.
E: Variable for fix langevin is invalid style
It must be an equal-style variable.
E: Fix langevin omega requires extended particles
One of the particles has radius 0.0.
E: Fix langevin angmom requires extended particles
This fix option cannot be used with point particles.
E: Cannot zero Langevin force of 0 atoms
The group has zero atoms, so you cannot request its force
be zeroed.
E: Fix langevin variable returned negative temperature
Self-explanatory.
E: Could not find fix_modify temperature ID
The compute ID for computing temperature does not exist.
E: Fix_modify temperature ID does not compute temperature
The compute ID assigned to the fix must compute temperature.
E: Fix langevin gjf cannot have period equal to dt/2
If the period is equal to dt/2 then division by zero can happen.
E: Fix langevin gjf should come before fix nve
Self-explanatory
E: Fix langevin gjf and respa are not compatible
Self-explanatory
W: Group for fix_modify temp != fix group
The fix_modify command is specifying a temperature computation that
computes a temperature on a different group of atoms than the fix
itself operates on. This is probably not what you want to do.
*/