Simplified GJF formalism (52a51ea4) · Commits · 郑智淋 / lammps

src/fix_langevin.cpp

+321 −443

Original line number	Diff line number	Diff line
		@@ -2,20 +2,16 @@
		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 authors: Carolyn Phillips (U Mich), reservoir energy tally
		Aidan Thompson (SNL) GJF formulation
		Charles Sievers (UC Davis) GJF-2GJ Implementation
		Niels Gronbech-Jensen (UC Davis) GJF-2GJ Formulation
		------------------------------------------------------------------------- */

		#include "fix_langevin.h"
		@@ -52,8 +48,7 @@ enum{CONSTANT,EQUAL,ATOM};
		FixLangevin::FixLangevin(LAMMPS lmp, int narg, char *arg) :
		Fix(lmp, narg, arg),
		gjfflag(0), gfactor1(NULL), gfactor2(NULL), ratio(NULL), tstr(NULL),
		flangevin(NULL), tforce(NULL), franprev(NULL), id_temp(NULL), random(NULL),
		lv(NULL), wildcard(NULL), bias(NULL)
		flangevin(NULL), tforce(NULL), franprev(NULL), id_temp(NULL), random(NULL), lv(NULL)
		{
		if (narg < 7) error->all(FLERR,"Illegal fix langevin command");

		@@ -98,7 +93,6 @@ FixLangevin::FixLangevin(LAMMPS lmp, int narg, char *arg) :
		oflag = 0;
		tallyflag = 0;
		zeroflag = 0;
		hsflag = 0;

		int iarg = 7;
		while (iarg < narg) {
		@@ -110,10 +104,7 @@ FixLangevin::FixLangevin(LAMMPS lmp, int narg, char *arg) :
		} else if (strcmp(arg[iarg],"gjf") == 0) {
		if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
		if (strcmp(arg[iarg+1],"no") == 0) gjfflag = 0;
		else if (strcmp(arg[iarg+1],"yes") == 0)
		error->all(FLERR,"GJF yes keyword is deprecated.\nPlease use vhalf or vfull.");
		else if (strcmp(arg[iarg+1],"vfull") == 0) {gjfflag = 1; hsflag = 0;}
		else if (strcmp(arg[iarg+1],"vhalf") == 0) {gjfflag = 1; hsflag = 1;}
		else if (strcmp(arg[iarg+1],"yes") == 0) gjfflag = 1;
		else error->all(FLERR,"Illegal fix langevin command");
		iarg += 2;
		} else if (strcmp(arg[iarg],"omega") == 0) {
		@@ -159,9 +150,7 @@ FixLangevin::FixLangevin(LAMMPS lmp, int narg, char *arg) :
		flangevin = NULL;
		flangevin_allocated = 0;
		franprev = NULL;
		wildcard = NULL;
		lv = NULL;
		bias = NULL;
		tforce = NULL;
		maxatom1 = maxatom2 = 0;

		@@ -170,8 +159,6 @@ FixLangevin::FixLangevin(LAMMPS lmp, int narg, char *arg) :
		// no need to set peratom_flag, b/c data is for internal use only

		if (gjfflag) {

		nvalues = 3;
		grow_arrays(atom->nmax);
		atom->add_callback(0);

		@@ -182,16 +169,11 @@ FixLangevin::FixLangevin(LAMMPS lmp, int narg, char *arg) :
		franprev[i][0] = 0.0;
		franprev[i][1] = 0.0;
		franprev[i][2] = 0.0;
		wildcard[i][0] = 0.0;
		wildcard[i][1] = 0.0;
		wildcard[i][2] = 0.0;
		if (hsflag) {
		lv[i][0] = 0.0;
		lv[i][1] = 0.0;
		lv[i][2] = 0.0;
		}
		}
		}

		}

		@@ -210,9 +192,7 @@ FixLangevin::~FixLangevin()

		if (gjfflag) {
		memory->destroy(franprev);
		memory->destroy(wildcard);
		if (hsflag) memory->destroy(lv);
		if (temperature && temperature->tempbias) memory->destroy(bias);
		memory->destroy(lv);
		atom->delete_callback(id,0);
		}
		}
		@@ -222,7 +202,8 @@ FixLangevin::~FixLangevin()
		int FixLangevin::setmask()
		{
		int mask = 0;
		if (gjfflag) mask \|= POST_INTEGRATE;
		if (gjfflag) mask \|= INITIAL_INTEGRATE;
		if (gjfflag) mask \|= INITIAL_INTEGRATE_RESPA;
		mask \|= POST_FORCE;
		mask \|= POST_FORCE_RESPA;
		mask \|= END_OF_STEP;
		@@ -234,14 +215,26 @@ int FixLangevin::setmask()

		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");

		// warn if any integrate fix comes after this one
		int before = 1;
		int flag = 0;
		for (int i = 0; i < modify->nfix; i++) {
		if (strcmp(id,modify->fix[i]->id) == 0) before = 0;
		else if ((modify->fmask[i] && strcmp(modify->fix[i]->style,"nve")==0) && before) flag = 1;
		}
		if (flag && comm->me == 0)
		error->all(FLERR,"Fix langevin gjf should come before fix nve");
		}

		if (oflag && !atom->sphere_flag)
		error->all(FLERR,"Fix langevin omega requires atom style sphere");
		if (ascale && !atom->ellipsoid_flag)
		error->all(FLERR,"Fix langevin angmom requires atom style ellipsoid");

		if (gjfflag && zeroflag && tallyflag)
		error->warning(FLERR,
		"Fix langevin: gjf, zero, and tally were all set correct energy tallying is not guaranteed");
		// check variable

		if (tstr) {
		@@ -281,16 +274,18 @@ void FixLangevin::init()
		error->one(FLERR,"Fix langevin angmom requires extended particles");
		}

		// set force prefactors

		if (!atom->rmass) {
		for (int i = 1; i <= atom->ntypes; i++) {
		gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
		if (gjfflag)
		if (!gjfflag)
		gfactor2[i] = sqrt(atom->mass[i]) *
		sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
		sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
		force->ftm2v;
		else
		gfactor2[i] = sqrt(atom->mass[i]) *
		sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
		sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
		force->ftm2v;
		gfactor1[i] *= 1.0/ratio[i];
		gfactor2[i] *= 1.0/sqrt(ratio[i]);
		@@ -303,161 +298,124 @@ void FixLangevin::init()
		if (strstr(update->integrate_style,"respa"))
		nlevels_respa = ((Respa *) update->integrate)->nlevels;

		if (strstr(update->integrate_style,"respa"))
		error->one(FLERR,"Fix langevin gjf not implemented with respa capabilities");

		if (gjfflag) gjffac = 1.0/sqrt(1.0+update->dt/2.0/t_period);

		if (gjfflag) gjfa = (1.0-update->dt/2.0/t_period)/(1.0+update->dt/2.0/t_period);
		if (gjfflag) gjfsib = sqrt(1.0+update->dt/2.0/t_period);
		}

		/* ---------------------------------------------------------------------- */

		void FixLangevin::setup(int vflag)
		{
		if (strstr(update->integrate_style,"verlet"))
		post_force(vflag);
		else {
		((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
		post_force_respa(vflag,nlevels_respa-1,0);
		((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
		}
		if (gjfflag){

		// update v of atoms in group
		double dtfm;
		double dt = update->dt;
		double **v = atom->v;
		double **f = atom->f;
		int *mask = atom->mask;
		int nlocal = atom->nlocal;
		if (igroup == atom->firstgroup) nlocal = atom->nfirst;
		double b[3] = {0.0,0.0,0.0};
		double *rmass = atom->rmass;
		double *mass = atom->mass;
		int *type = atom->type;
		if (rmass) {
		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit) {
		dtfm = 0.5 * dt / rmass[i];
		v[i][0] -= dtfm * f[i][0];
		v[i][1] -= dtfm * f[i][1];
		v[i][2] -= dtfm * f[i][2];
		if (tbiasflag)
		temperature->remove_bias(i,v[i]);
		v[i][0] /= gjfagjfsibgjfsib;
		v[i][1] /= gjfagjfsibgjfsib;
		v[i][2] /= gjfagjfsibgjfsib;
		if (tbiasflag)
		temperature->restore_bias(i,v[i]);
		}

		for (int i = 0; i < nlocal; i++) {
		f[i][0] = wildcard[i][0];
		f[i][1] = wildcard[i][1];
		f[i][2] = wildcard[i][2];
		b[0] = v[i][0];
		b[1] = v[i][1];
		b[2] = v[i][2];
		if (tbiasflag == BIAS) temperature->remove_bias(i,v[i]);
		wildcard[i][0] = v[i][0];
		wildcard[i][1] = v[i][1];
		wildcard[i][2] = v[i][2];
		if (tbiasflag == BIAS) {
		} else {
		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit) {
		dtfm = 0.5 * dt / mass[type[i]];
		v[i][0] -= dtfm * f[i][0];
		v[i][1] -= dtfm * f[i][1];
		v[i][2] -= dtfm * f[i][2];
		if (tbiasflag)
		temperature->remove_bias(i,v[i]);
		v[i][0] /= gjfagjfsibgjfsib;
		v[i][1] /= gjfagjfsibgjfsib;
		v[i][2] /= gjfagjfsibgjfsib;
		if (tbiasflag)
		temperature->restore_bias(i,v[i]);
		bias[i][0] = b[0] - wildcard[i][0];
		bias[i][1] = b[1] - wildcard[i][1];
		bias[i][2] = b[2] - wildcard[i][2];
		}
		}
		}
		if (strstr(update->integrate_style,"verlet"))
		post_force(vflag);
		else {
		((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
		post_force_respa(vflag,nlevels_respa-1,0);
		((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
		}

		/* ----------------------------------------------------------------------
		allow for both per-type and per-atom mass
		------------------------------------------------------------------------- */

		void FixLangevin::post_integrate()
		{
		if (gjfflag){
		double dtfm;
		double dt = update->dt;
		double dtf = 0.5 * dt * force->ftm2v;

		// update v of atoms in group

		double **x = atom->x;
		double **v = atom->v;
		double **f = atom->f;
		double **v = atom->v;
		int *mask = atom->mask;
		int nlocal = atom->nlocal;
		double *rmass = atom->rmass;
		double *mass = atom->mass;
		int *type = atom->type;
		int *mask = atom->mask;
		int nlocal = atom->nlocal;
		if (igroup == atom->firstgroup) nlocal = atom->nfirst;

		// zero option
		double vsum[3],vsumall[3];
		bigint count;

		if (zeroflag) {
		vsum[0] = vsum[1] = vsum[2] = 0.0;
		count = group->count(igroup);
		if (count == 0)
		error->all(FLERR,"Cannot zero Langevin force of 0 atoms");
		}

		if (rmass) {
		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit) {
		dtfm = dtf / rmass[i];
		x[i][0] -= dt * v[i][0];
		x[i][1] -= dt * v[i][1];
		x[i][2] -= dt * v[i][2];
		v[i][0] = gjffac * (wildcard[i][0] + dtfm * franprev[i][0] + dtfm * f[i][0]);
		v[i][1] = gjffac * (wildcard[i][1] + dtfm * franprev[i][1] + dtfm * f[i][1]);
		v[i][2] = gjffac * (wildcard[i][2] + dtfm * franprev[i][2] + dtfm * f[i][2]);
		if (tbiasflag == BIAS)
		for (int j = 0; j < 3; j++) {
		if (wildcard[i][j] == 0) {
		v[i][j] /= gjffac * gjffac;
		}
		}
		x[i][0] += gjffac * dt * v[i][0];
		x[i][1] += gjffac * dt * v[i][1];
		x[i][2] += gjffac * dt * v[i][2];
		if (tbiasflag == BIAS)
		for (int j = 0; j < 3; j++) {
		if (wildcard[i][j] == 0)
		v[i][j] *= gjffac;
		v[i][j] += bias[i][j];
		x[i][j] += dt * bias[i][j];
		}
		}

		dtfm = 0.5 * dt / rmass[i];
		v[i][0] += dtfm * f[i][0];
		v[i][1] += dtfm * f[i][1];
		v[i][2] += dtfm * f[i][2];
		lv[i][0] = f[i][0];
		lv[i][1] = f[i][1];
		lv[i][2] = f[i][2];
		}
		//
		} else {
		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit) {
		dtfm = dtf / mass[type[i]];
		x[i][0] -= dt * v[i][0];
		x[i][1] -= dt * v[i][1];
		x[i][2] -= dt * v[i][2];
		v[i][0] = gjffac * (wildcard[i][0] + dtfm * franprev[i][0] + dtfm * f[i][0]);
		v[i][1] = gjffac * (wildcard[i][1] + dtfm * franprev[i][1] + dtfm * f[i][1]);
		v[i][2] = gjffac * (wildcard[i][2] + dtfm * franprev[i][2] + dtfm * f[i][2]);
		if (tbiasflag == BIAS)
		for (int j = 0; j < 3; j++) {
		if (wildcard[i][j] == 0) {
		v[i][j] /= gjffac*gjffac;
		}
		dtfm = 0.5 * dt / mass[type[i]];
		v[i][0] += dtfm * f[i][0];
		v[i][1] += dtfm * f[i][1];
		v[i][2] += dtfm * f[i][2];
		lv[i][0] = v[i][0];
		lv[i][1] = v[i][1];
		lv[i][2] = v[i][2];
		}
		x[i][0] += gjffac * dt * v[i][0];
		x[i][1] += gjffac * dt * v[i][1];
		x[i][2] += gjffac * dt * v[i][2];
		if (tbiasflag == BIAS)
		for (int j = 0; j < 3; j++) {
		if (wildcard[i][j] == 0)
		v[i][j] *= gjffac;
		v[i][j] += bias[i][j];
		x[i][j] += dt * bias[i][j];
		}
		if (zeroflag){
		vsum[0] += gjffac * dtfm * franprev[i][0];
		vsum[1] += gjffac * dtfm * franprev[i][1];
		vsum[2] += gjffac * dtfm * franprev[i][2];
		}
		}

		/* ---------------------------------------------------------------------- */

		void FixLangevin::initial_integrate_respa(int vflag, int ilevel, int /* iloop */){
		if (ilevel == respa_level-1) initial_integrate(vflag);
		}

		if (zeroflag) {
		MPI_Allreduce(vsum,vsumall,3,MPI_DOUBLE,MPI_SUM,world);
		vsumall[0] /= count;
		vsumall[1] /= count;
		vsumall[2] /= count;
		for (int i = 0; i < nlocal; i++) {
		/* ---------------------------------------------------------------------- */

		void FixLangevin::initial_integrate(int /* vflag */)
		{
		double **v = atom->v;
		double **f = atom->f;
		int *mask = atom->mask;
		int nlocal = atom->nlocal;

		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit){
		v[i][0] -= vsumall[0];
		v[i][1] -= vsumall[1];
		v[i][2] -= vsumall[2];
		}
		}
		f[i][0] /= gjfa;
		f[i][1] /= gjfa;
		f[i][2] /= gjfa;
		v[i][0] = lv[i][0];
		v[i][1] = lv[i][1];
		v[i][2] = lv[i][2];
		}
		}

		@@ -644,8 +602,9 @@ void FixLangevin::post_force_templated()
		// sum random force over all atoms in group
		// subtract sum/count from each atom in group

		double fdrag[3],fran[3],fsum[3],fsumall[3], rantemp[3];
		double fdrag[3],fran[3],fsum[3],fsumall[3];
		bigint count;
		double fswap;

		double boltz = force->boltz;
		double dt = update->dt;
		@@ -672,18 +631,17 @@ void FixLangevin::post_force_templated()
		flangevin_allocated = 1;
		}

		if (Tp_BIAS && !gjfflag) temperature->compute_scalar();
		else if (Tp_BIAS && update->ntimestep == update->beginstep && gjfflag) temperature->compute_scalar();
		if (Tp_BIAS) temperature->compute_scalar();

		for (int i = 0; i < nlocal; i++) {
		if (mask[i] & groupbit) {
		if (Tp_TSTYLEATOM) tsqrt = sqrt(tforce[i]);
		if (Tp_RMASS) {
		gamma1 = -rmass[i] / t_period / ftm2v;
		if (Tp_GJF)
		gamma2 = sqrt(rmass[i]) * sqrt(2.0 * boltz / t_period / dt / mvv2e) / ftm2v;
		else
		if (!Tp_GJF)
		gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
		else
		gamma2 = sqrt(rmass[i]) * sqrt(2.0*boltz/t_period/dt/mvv2e) / ftm2v;
		gamma1 *= 1.0/ratio[type[i]];
		gamma2 = 1.0/sqrt(ratio[type[i]]) tsqrt;
		} else {
		@@ -691,15 +649,16 @@ void FixLangevin::post_force_templated()
		gamma2 = gfactor2[type[i]] * tsqrt;
		}

		if (Tp_GJF) {
		fran[0] = gamma2 * random->gaussian();
		fran[1] = gamma2 * random->gaussian();
		fran[2] = gamma2 * random->gaussian();
		} else {
		if (!Tp_GJF){
		fran[0] = gamma2*(random->uniform()-0.5);
		fran[1] = gamma2*(random->uniform()-0.5);
		fran[2] = gamma2*(random->uniform()-0.5);
		}
		else{
		fran[0] = gamma2*random->gaussian();
		fran[1] = gamma2*random->gaussian();
		fran[2] = gamma2*random->gaussian();
		}

		if (Tp_BIAS) {
		temperature->remove_bias(i,v[i]);
		@@ -717,21 +676,35 @@ void FixLangevin::post_force_templated()
		}

		if (Tp_GJF) {
		wildcard[i][0] = f[i][0];
		wildcard[i][1] = f[i][1];
		wildcard[i][2] = f[i][2];

		rantemp[0] = fran[0];
		rantemp[1] = fran[1];
		rantemp[2] = fran[2];

		fran[0] = franprev[i][0];
		fran[1] = franprev[i][1];
		fran[2] = franprev[i][2];

		fdrag[0] = -2t_period((2gjffac)-(1.0/gjffac)-1.0)/dt;
		fdrag[1] = -2t_period((2gjffac)-(1.0/gjffac)-1.0)/dt;
		fdrag[2] = -2t_period((2gjffac)-(1.0/gjffac)-1.0)/dt;
		if (Tp_BIAS)
		temperature->remove_bias(i,v[i]);
		lv[i][0] = gjfsib*v[i][0];
		lv[i][1] = gjfsib*v[i][1];
		lv[i][2] = gjfsib*v[i][2];
		if (Tp_BIAS)
		temperature->restore_bias(i,v[i]);
		if (Tp_BIAS)
		temperature->restore_bias(i,lv[i]);

		fswap = 0.5*(fran[0]+franprev[i][0]);
		franprev[i][0] = fran[0];
		fran[0] = fswap;
		fswap = 0.5*(fran[1]+franprev[i][1]);
		franprev[i][1] = fran[1];
		fran[1] = fswap;
		fswap = 0.5*(fran[2]+franprev[i][2]);
		franprev[i][2] = fran[2];
		fran[2] = fswap;

		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];
		@@ -739,62 +712,16 @@ void FixLangevin::post_force_templated()
		f[i][2] += fdrag[2] + fran[2];

		if (Tp_TALLY) {
		if (Tp_GJF && update->ntimestep != update->beginstep){
		if (Tp_BIAS) {
		temperature->remove_bias(i,v[i]);
		fdrag[0] = gamma1gjffacgjffac*v[i][0];
		fdrag[1] = gamma1gjffacgjffac*v[i][1];
		fdrag[2] = gamma1gjffacgjffac*v[i][2];
		temperature->restore_bias(i,v[i]);
		} else {
		fdrag[0] = gamma1gjffacgjffac*v[i][0];
		fdrag[1] = gamma1gjffacgjffac*v[i][1];
		fdrag[2] = gamma1gjffacgjffac*v[i][2];
		}
		fran[0] *= gjffac;
		fran[1] *= gjffac;
		fran[2] *= gjffac;
		}
		else if (Tp_GJF && update->ntimestep == update->beginstep){
		fdrag[0] = 0.0;
		fdrag[1] = 0.0;
		fdrag[2] = 0.0;
		}
		flangevin[i][0] = fdrag[0] + fran[0];
		flangevin[i][1] = fdrag[1] + fran[1];
		flangevin[i][2] = fdrag[2] + fran[2];
		}

		if (Tp_ZERO) {
		if (!Tp_GJF){
		fsum[0] += fran[0];
		fsum[1] += fran[1];
		fsum[2] += fran[2];
		}
		else {
		fsum[0] += franprev[i][0];
		fsum[1] += franprev[i][1];
		fsum[2] += franprev[i][2];
		}
		}

		if (Tp_GJF)
		{
		franprev[i][0] = rantemp[0];
		franprev[i][1] = rantemp[1];
		franprev[i][2] = rantemp[2];

		if (hsflag){
		lv[i][0] = v[i][0];
		lv[i][1] = v[i][1];
		lv[i][2] = v[i][2];
		if (tbiasflag == BIAS) {
		lv[i][0] += bias[i][0];
		lv[i][1] += bias[i][1];
		lv[i][2] += bias[i][2];
		}
		}
		}
		}
		}

		@@ -977,47 +904,35 @@ void FixLangevin::end_of_step()
		if (!tallyflag && !gjfflag) return;

		double **v = atom->v;
		double **f = atom->f;
		int *mask = atom->mask;
		int nlocal = atom->nlocal;
		double b[3] = {0.0,0.0,0.0};

		if (gjfflag && tbiasflag == BIAS) temperature->compute_scalar();

		energy_onestep = 0.0;

		if (gjfflag){
		double tmp[3];
		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit){
		if (gjfflag){
		b[0] = v[i][0];
		b[1] = v[i][1];
		b[2] = v[i][2];
		f[i][0] = wildcard[i][0];
		f[i][1] = wildcard[i][1];
		f[i][2] = wildcard[i][2];
		if (tbiasflag == BIAS) temperature->remove_bias(i,v[i]);
		wildcard[i][0] = v[i][0];
		wildcard[i][1] = v[i][1];
		wildcard[i][2] = v[i][2];
		if (tbiasflag == BIAS) {
		bias[i][0] = b[0] - v[i][0];
		bias[i][1] = b[1] - v[i][1];
		bias[i][2] = b[2] - v[i][2];
		temperature->restore_bias(i, v[i]);
		}
		if (hsflag){
		tmp[0] = v[i][0];
		tmp[1] = v[i][1];
		tmp[2] = v[i][2];
		v[i][0] = lv[i][0];
		v[i][1] = lv[i][1];
		v[i][2] = lv[i][2];
		lv[i][0] = tmp[0];
		lv[i][1] = tmp[1];
		lv[i][2] = tmp[2];
		}
		}

		if (tallyflag)
		for (int i = 0; i < nlocal; i++)
		if (mask[i] & groupbit)
		energy_onestep += flangevin[i][0]v[i][0] + flangevin[i][1]v[i][1] +
		flangevin[i][2]*v[i][2];
		}
		if (tallyflag) {

		energy += energy_onestep*update->dt;
		}
		}

		/* ---------------------------------------------------------------------- */

		@@ -1070,7 +985,7 @@ int FixLangevin::modify_param(int narg, char **arg)

		double FixLangevin::compute_scalar()
		{
		if (!tallyflag && !flangevin_allocated) return 0.0;
		if (!tallyflag \|\| !flangevin_allocated) return 0.0;

		// capture the very first energy transfer to thermal reservoir

		@@ -1088,11 +1003,8 @@ double FixLangevin::compute_scalar()
		}

		// convert midstep energy back to previous fullstep energy
		double energy_me;
		if (gjfflag)
		energy_me = energy - energy_onestep*update->dt;
		else
		energy_me = energy - 0.5energy_onestepupdate->dt;

		double energy_me = energy - 0.5energy_onestepupdate->dt;

		double energy_all;
		MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
		@@ -1119,9 +1031,7 @@ void FixLangevin::extract(const char str, int &dim)
		double FixLangevin::memory_usage()
		{
		double bytes = 0.0;
		if (gjfflag) bytes += atom->nmax32 * sizeof(double);
		if (gjfflag) if (hsflag) bytes += atom->nmax3 sizeof(double);
		if (gjfflag && tbiasflag == BIAS) bytes += atom->nmax3 sizeof(double);
		if (gjfflag) bytes += atom->nmax6 sizeof(double);
		if (tallyflag) bytes += atom->nmax3 sizeof(double);
		if (tforce) bytes += atom->nmax * sizeof(double);
		return bytes;
		@@ -1134,9 +1044,7 @@ double FixLangevin::memory_usage()
		void FixLangevin::grow_arrays(int nmax)
		{
		memory->grow(franprev,nmax,3,"fix_langevin:franprev");
		memory->grow(wildcard,nmax,3,"fix_langevin:wildcard");
		if (hsflag) memory->grow(lv,nmax,3,"fix_langevin:lv");
		if (tbiasflag == BIAS) memory->grow(bias,nmax,3,"fix_langevin:bias");
		memory->grow(lv,nmax,3,"fix_langevin:lv");
		}

		/* ----------------------------------------------------------------------
		@@ -1148,20 +1056,10 @@ void FixLangevin::copy_arrays(int i, int j, int /delflag/)
		franprev[j][0] = franprev[i][0];
		franprev[j][1] = franprev[i][1];
		franprev[j][2] = franprev[i][2];
		wildcard[j][0] = wildcard[i][0];
		wildcard[j][1] = wildcard[i][1];
		wildcard[j][2] = wildcard[i][2];
		if (hsflag) {
		lv[j][0] = lv[i][0];
		lv[j][1] = lv[i][1];
		lv[j][2] = lv[i][2];
		}
		if (tbiasflag == BIAS){
		bias[j][0] = bias[i][0];
		bias[j][1] = bias[i][1];
		bias[j][2] = bias[i][2];
		}
		}

		/* ----------------------------------------------------------------------
		pack values in local atom-based array for exchange with another proc
		@@ -1173,19 +1071,9 @@ int FixLangevin::pack_exchange(int i, double *buf)
		buf[n++] = franprev[i][0];
		buf[n++] = franprev[i][1];
		buf[n++] = franprev[i][2];
		buf[n++] = wildcard[i][0];
		buf[n++] = wildcard[i][1];
		buf[n++] = wildcard[i][2];
		if (hsflag){
		buf[n++] = lv[i][0];
		buf[n++] = lv[i][1];
		buf[n++] = lv[i][2];
		}
		if (tbiasflag == BIAS){
		buf[n++] = bias[i][0];
		buf[n++] = bias[i][1];
		buf[n++] = bias[i][2];
		}
		return n;
		}

		@@ -1199,18 +1087,8 @@ int FixLangevin::unpack_exchange(int nlocal, double *buf)
		franprev[nlocal][0] = buf[n++];
		franprev[nlocal][1] = buf[n++];
		franprev[nlocal][2] = buf[n++];
		wildcard[nlocal][0] = buf[n++];
		wildcard[nlocal][1] = buf[n++];
		wildcard[nlocal][2] = buf[n++];
		if (hsflag){
		lv[nlocal][0] = buf[n++];
		lv[nlocal][1] = buf[n++];
		lv[nlocal][2] = buf[n++];
		}
		if (tbiasflag == BIAS){
		bias[nlocal][0] = buf[n++];
		bias[nlocal][1] = buf[n++];
		bias[nlocal][2] = buf[n++];
		}
		return n;
		}
		No newline at end of file

src/fix_langevin.h

+58 −91

File changed.

Preview size limit exceeded, changes collapsed.

Admin message