Merge pull request #806 from vishalkenchan/new_pair_lj_cut_coul_wolf (9cf99ef4) · Commits · 郑智淋 / lammps

doc/src/pair_lj.txt

+28 −0

Original line number	Diff line number	Diff line
		@@ -32,6 +32,8 @@ pair_style lj/cut/coul/long/omp command :h3
		pair_style lj/cut/coul/msm command :h3
		pair_style lj/cut/coul/msm/gpu command :h3
		pair_style lj/cut/coul/msm/omp command :h3
		pair_style lj/cut/coul/wolf command :h3
		pair_style lj/cut/coul/wolf/omp command :h3
		pair_style lj/cut/tip4p/cut command :h3
		pair_style lj/cut/tip4p/cut/omp command :h3
		pair_style lj/cut/tip4p/long command :h3
		@@ -63,6 +65,10 @@ args = list of arguments for a particular style :ul
		{lj/cut/coul/msm} args = cutoff (cutoff2)
		cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
		cutoff2 = global cutoff for Coulombic (optional) (distance units)
		{lj/cut/coul/wolf} args = alpha cutoff (cutoff2)
		alpha = damping parameter (inverse distance units)
		cutoff = global cutoff for LJ (and Coulombic if only 2 arg) (distance units)
		cutoff2 = global cutoff for Coulombic (optional) (distance units)
		{lj/cut/tip4p/cut} args = otype htype btype atype qdist cutoff (cutoff2)
		otype,htype = atom types for TIP4P O and H
		btype,atype = bond and angle types for TIP4P waters
		@@ -115,6 +121,10 @@ pair_style lj/cut/tip4p/cut 1 2 7 8 0.15 12.0 10.0
		pair_coeff * * 100.0 3.0
		pair_coeff 1 1 100.0 3.5 9.0 :pre

		pair_style lj/cut/coul/wolf 0.2 5. 10.0
		pair_coeff * * 1.0 1.0
		pair_coeff 1 1 1.0 1.0 2.5 :pre

		pair_style lj/cut/tip4p/long 1 2 7 8 0.15 12.0
		pair_style lj/cut/tip4p/long 1 2 7 8 0.15 12.0 10.0
		pair_coeff * * 100.0 3.0
		@@ -179,6 +189,24 @@ that a term is added for the "core/shell
		model"_Section_howto.html#howto_25 to allow charges on core and shell
		particles to be separated by r = 0.0.

		Style {coul/wolf} adds a Coulombic pairwise interaction via the Wolf
		summation method, described in "Wolf"_#Wolf1, given by:

		:c,image(Eqs/pair_coul_wolf.jpg)

		where {alpha} is the damping parameter, and erfc() is the
		complementary error-function terms. This potential
		is essentially a short-range, spherically-truncated,
		charge-neutralized, shifted, pairwise {1/r} summation. With a
		manipulation of adding and subtracting a self term (for i = j) to the
		first and second term on the right-hand-side, respectively, and a
		small enough {alpha} damping parameter, the second term shrinks and
		the potential becomes a rapidly-converging real-space summation. With
		a long enough cutoff and small enough alpha parameter, the energy and
		forces calculated by the Wolf summation method approach those of the
		Ewald sum. So it is a means of getting effective long-range
		interactions with a short-range potential.

		Styles {lj/cut/tip4p/cut} and {lj/cut/tip4p/long} implement the TIP4P
		water model of "(Jorgensen)"_#Jorgensen2, which introduces a massless
		site located a short distance away from the oxygen atom along the

src/USER-OMP/pair_lj_cut_coul_wolf_omp.cpp

0 → 100644

+204 −0

Original line number	Diff line number	Diff line
		/* ----------------------------------------------------------------------
		LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
		http://lammps.sandia.gov, Sandia National Laboratories
		Steve Plimpton, sjplimp@sandia.gov

		This software is distributed under the GNU General Public License.

		See the README file in the top-level LAMMPS directory.
		------------------------------------------------------------------------- */

		/* ----------------------------------------------------------------------
		Contributing author: Axel Kohlmeyer (Temple U)
		------------------------------------------------------------------------- */

		#include <math.h>
		#include "pair_lj_cut_coul_wolf_omp.h"
		#include "atom.h"
		#include "comm.h"
		#include "force.h"
		#include "neighbor.h"
		#include "neigh_list.h"

		#include "suffix.h"
		#include "math_const.h"
		using namespace LAMMPS_NS;
		using namespace MathConst;

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

		PairLJCutCoulWolfOMP::PairLJCutCoulWolfOMP(LAMMPS *lmp) :
		PairLJCutCoulWolf(lmp), ThrOMP(lmp, THR_PAIR)
		{
		suffix_flag \|= Suffix::OMP;
		respa_enable = 0;
		}

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

		void PairLJCutCoulWolfOMP::compute(int eflag, int vflag)
		{
		if (eflag \|\| vflag) {
		ev_setup(eflag,vflag);
		} else evflag = vflag_fdotr = 0;

		const int nall = atom->nlocal + atom->nghost;
		const int nthreads = comm->nthreads;
		const int inum = list->inum;

		#if defined(_OPENMP)
		#pragma omp parallel default(none) shared(eflag,vflag)
		#endif
		{
		int ifrom, ito, tid;

		loop_setup_thr(ifrom, ito, tid, inum, nthreads);
		ThrData *thr = fix->get_thr(tid);
		thr->timer(Timer::START);
		ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);

		if (evflag) {
		if (eflag) {
		if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
		else eval<1,1,0>(ifrom, ito, thr);
		} else {
		if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
		else eval<1,0,0>(ifrom, ito, thr);
		}
		} else {
		if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
		else eval<0,0,0>(ifrom, ito, thr);
		}

		thr->timer(Timer::PAIR);
		reduce_thr(this, eflag, vflag, thr);
		} // end of omp parallel region
		}

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

		template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
		void PairLJCutCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
		{
		int i,j,ii,jj,jnum,itype,jtype;
		double qitmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
		double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
		double prefactor,erfcc,erfcd,v_sh,dvdrr,e_self,qisq;
		int ilist,jlist,numneigh,*firstneigh;

		evdwl = ecoul = 0.0;

		const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
		dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
		const double * _noalias const q = atom->q;
		const int * _noalias const type = atom->type;
		const int nlocal = atom->nlocal;
		const double * _noalias const special_coul = force->special_coul;
		const double * _noalias const special_lj = force->special_lj;
		const double qqrd2e = force->qqrd2e;
		double fxtmp,fytmp,fztmp;

		// self and shifted Coulombic energy

		e_self = v_sh = 0.0;
		const double e_shift = erfc(alf*cut_coul)/cut_coul;
		const double f_shift = -(e_shift+ 2.0*alf/MY_PIS
		* exp(-alfalfcut_coul*cut_coul)) / cut_coul;

		ilist = list->ilist;
		numneigh = list->numneigh;
		firstneigh = list->firstneigh;

		// loop over neighbors of my atoms

		for (ii = iifrom; ii < iito; ++ii) {

		i = ilist[ii];
		qitmp = q[i];
		xtmp = x[i].x;
		ytmp = x[i].y;
		ztmp = x[i].z;
		itype = type[i];
		jlist = firstneigh[i];
		jnum = numneigh[i];
		fxtmp=fytmp=fztmp=0.0;

		if (EFLAG) {
		e_self = -(e_shift/2.0 + alf/MY_PIS) * qitmpqitmpqqrd2e;
		ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr);
		}

		for (jj = 0; jj < jnum; jj++) {
		j = jlist[jj];
		factor_lj = special_lj[sbmask(j)];
		factor_coul = special_coul[sbmask(j)];
		j &= NEIGHMASK;

		delx = xtmp - x[j].x;
		dely = ytmp - x[j].y;
		delz = ztmp - x[j].z;
		rsq = delxdelx + delydely + delz*delz;
		jtype = type[j];

		if (rsq < cutsq[itype][jtype]) {
		r2inv = 1.0/rsq;

		if (rsq < cut_coulsq) {
		r = sqrt(rsq);
		prefactor = qqrd2eqitmpq[j]/r;
		erfcc = erfc(alf*r);
		erfcd = exp(-alfalfr*r);
		v_sh = (erfcc - e_shiftr) prefactor;
		dvdrr = (erfcc/rsq + 2.0alf/MY_PIS erfcd/r) + f_shift;
		forcecoul = dvdrrrsqprefactor;
		if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
		} else forcecoul = 0.0;

		if (rsq < cut_ljsq[itype][jtype]) {
		r6inv = r2invr2invr2inv;
		forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
		} else forcelj = 0.0;

		fpair = (forcecoul + factor_ljforcelj) r2inv;

		fxtmp += delx*fpair;
		fytmp += dely*fpair;
		fztmp += delz*fpair;
		if (NEWTON_PAIR \|\| j < nlocal) {
		f[j].x -= delx*fpair;
		f[j].y -= dely*fpair;
		f[j].z -= delz*fpair;
		}

		if (EFLAG) {
		if (rsq < cut_ljsq[itype][jtype]) {
		evdwl = r6inv(lj3[itype][jtype]r6inv-lj4[itype][jtype]) -
		offset[itype][jtype];
		evdwl *= factor_lj;
		} else evdwl = 0.0;

		if (rsq < cut_coulsq) {
		ecoul = v_sh;
		if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
		} else ecoul = 0.0;
		}

		if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
		evdwl,ecoul,fpair,delx,dely,delz,thr);
		}
		}
		f[i].x += fxtmp;
		f[i].y += fytmp;
		f[i].z += fztmp;
		}
		}

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

		double PairLJCutCoulWolfOMP::memory_usage()
		{
		double bytes = memory_usage_thr();
		bytes += PairLJCutCoulWolf::memory_usage();

		return bytes;
		}

src/USER-OMP/pair_lj_cut_coul_wolf_omp.h

0 → 100644

+48 −0

Original line number	Diff line number	Diff line
		/* -- c++ -- ----------------------------------------------------------
		LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
		http://lammps.sandia.gov, Sandia National Laboratories
		Steve Plimpton, sjplimp@sandia.gov

		Copyright (2003) Sandia Corporation. Under the terms of Contract
		DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
		certain rights in this software. This software is distributed under
		the GNU General Public License.

		See the README file in the top-level LAMMPS directory.
		------------------------------------------------------------------------- */

		/* ----------------------------------------------------------------------
		Contributing author: Axel Kohlmeyer (Temple U)
		------------------------------------------------------------------------- */

		#ifdef PAIR_CLASS

		PairStyle(lj/cut/coul/wolf/omp,PairLJCutCoulWolfOMP)

		#else

		#ifndef LMP_PAIR_LJ_CUT_COUL_WOLF_OMP_H
		#define LMP_PAIR_LJ_CUT_COUL_WOLF_OMP_H

		#include "pair_lj_cut_coul_wolf.h"
		#include "thr_omp.h"

		namespace LAMMPS_NS {

		class PairLJCutCoulWolfOMP : public PairLJCutCoulWolf, public ThrOMP {

		public:
		PairLJCutCoulWolfOMP(class LAMMPS *);

		virtual void compute(int, int);
		virtual double memory_usage();

		private:
		template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
		void eval(int ifrom, int ito, ThrData * const thr);
		};

		}

		#endif
		#endif

src/pair_lj_cut_coul_wolf.cpp

0 → 100644

+447 −0

Original line number	Diff line number	Diff line
		/* ----------------------------------------------------------------------
		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: Vishal Boddu (FAU)
		------------------------------------------------------------------------- */

		#include <math.h>
		#include <stdio.h>
		#include <stdlib.h>
		#include <string.h>
		#include "pair_lj_cut_coul_wolf.h"
		#include "atom.h"
		#include "comm.h"
		#include "force.h"
		#include "neighbor.h"
		#include "neigh_list.h"
		#include "neigh_request.h"
		#include "update.h"
		#include "integrate.h"
		#include "math_const.h"
		#include "memory.h"
		#include "error.h"

		using namespace LAMMPS_NS;
		using namespace MathConst;

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

		PairLJCutCoulWolf::PairLJCutCoulWolf(LAMMPS *lmp) : Pair(lmp)
		{
		single_enable = 0;
		writedata = 1;
		}

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

		PairLJCutCoulWolf::~PairLJCutCoulWolf()
		{
		if (allocated) {
		memory->destroy(setflag);
		memory->destroy(cutsq);

		memory->destroy(cut_lj);
		memory->destroy(cut_ljsq);
		memory->destroy(epsilon);
		memory->destroy(sigma);
		memory->destroy(lj1);
		memory->destroy(lj2);
		memory->destroy(lj3);
		memory->destroy(lj4);
		memory->destroy(offset);
		}
		}

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

		void PairLJCutCoulWolf::compute(int eflag, int vflag)
		{
		int i,j,ii,jj,inum,jnum,itype,jtype;
		double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
		double rsq,r2inv,r6inv,forcelj,factor_lj,forcecoul,factor_coul;
		double prefactor;
		double r;
		int ilist,jlist,numneigh,*firstneigh;
		double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq;

		evdwl = 0.0;
		ecoul = 0.0;
		if (eflag \|\| vflag) ev_setup(eflag,vflag);
		else evflag = vflag_fdotr = 0;

		double **x = atom->x;
		double **f = atom->f;
		double *q = atom->q;
		int *type = atom->type;
		int nlocal = atom->nlocal;
		double *special_lj = force->special_lj;
		double *special_coul = force->special_coul;
		int newton_pair = force->newton_pair;
		double qqrd2e = force->qqrd2e;

		// self and shifted Coulombic energy

		e_self = v_sh = 0.0;
		e_shift = erfc(alf*cut_coul)/cut_coul;
		f_shift = -(e_shift+ 2.0*alf/MY_PIS
		* exp(-alfalfcut_coul*cut_coul)) / cut_coul;

		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];
		qtmp = q[i];
		xtmp = x[i][0];
		ytmp = x[i][1];
		ztmp = x[i][2];
		itype = type[i];
		jlist = firstneigh[i];
		jnum = numneigh[i];

		qisq = qtmp*qtmp;
		e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e;
		if (eflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);

		for (jj = 0; jj < jnum; jj++) {
		j = jlist[jj];
		factor_lj = special_lj[sbmask(j)];
		factor_coul = special_coul[sbmask(j)];
		j &= NEIGHMASK;

		delx = xtmp - x[j][0];
		dely = ytmp - x[j][1];
		delz = ztmp - x[j][2];
		rsq = delxdelx + delydely + delz*delz;
		jtype = type[j];

		if (rsq < cutsq[itype][jtype]) {
		r2inv = 1.0/rsq;

		if (rsq < cut_coulsq) {
		r = sqrt(rsq);
		prefactor = qqrd2eqtmpq[j]/r;
		erfcc = erfc(alf*r);
		erfcd = exp(-alfalfr*r);
		v_sh = (erfcc - e_shiftr) prefactor;
		dvdrr = (erfcc/rsq + 2.0alf/MY_PIS erfcd/r) + f_shift;
		forcecoul = dvdrrrsqprefactor;
		if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
		} else forcecoul = 0.0;

		if (rsq < cut_ljsq[itype][jtype]) {
		r6inv = r2invr2invr2inv;
		forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
		} else forcelj = 0.0;

		fpair = (forcecoul + factor_ljforcelj) r2inv;

		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) {
		if (rsq < cut_ljsq[itype][jtype]) {
		evdwl = r6inv(lj3[itype][jtype]r6inv-lj4[itype][jtype]) -
		offset[itype][jtype];
		evdwl *= factor_lj;
		} else evdwl = 0.0;

		if (rsq < cut_coulsq) {
		ecoul = v_sh;
		if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
		} else ecoul = 0.0;
		}

		if (evflag) ev_tally(i,j,nlocal,newton_pair,
		evdwl,ecoul,fpair,delx,dely,delz);
		}
		}
		}

		if (vflag_fdotr) virial_fdotr_compute();
		}

		/* ----------------------------------------------------------------------
		allocate all arrays
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::allocate()
		{
		allocated = 1;
		int n = atom->ntypes;

		memory->create(setflag,n+1,n+1,"pair:setflag");
		for (int i = 1; i <= n; i++)
		for (int j = i; j <= n; j++)
		setflag[i][j] = 0;

		memory->create(cutsq,n+1,n+1,"pair:cutsq");
		memory->create(cut_lj,n+1,n+1,"pair:cut");
		memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
		memory->create(epsilon,n+1,n+1,"pair:epsilon");
		memory->create(sigma,n+1,n+1,"pair:sigma");
		memory->create(lj1,n+1,n+1,"pair:lj1");
		memory->create(lj2,n+1,n+1,"pair:lj2");
		memory->create(lj3,n+1,n+1,"pair:lj3");
		memory->create(lj4,n+1,n+1,"pair:lj4");
		memory->create(offset,n+1,n+1,"pair:offset");
		}

		/* ----------------------------------------------------------------------
		global settings
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::settings(int narg, char **arg)
		{
		if (narg < 2 \|\| narg > 3) error->all(FLERR,"Illegal pair_style command");

		alf = force->numeric(FLERR,arg[0]);
		cut_lj_global = force->numeric(FLERR,arg[1]);
		if (narg == 2) cut_coul = cut_lj_global;

		if (allocated) {
		int i,j;
		for (i = 1; i <= atom->ntypes; i++)
		for (j = i; j <= atom->ntypes; j++)
		if (setflag[i][j]) cut_lj[i][j] = cut_lj_global;
		}
		}

		/* ----------------------------------------------------------------------
		set coeffs for one or more type pairs
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::coeff(int narg, char **arg)
		{
		if (narg < 4 \|\| narg > 5)
		error->all(FLERR,"Incorrect args for pair coefficients");
		if (!allocated) allocate();

		int ilo,ihi,jlo,jhi;
		force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
		force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);

		double epsilon_one = force->numeric(FLERR,arg[2]);
		double sigma_one = force->numeric(FLERR,arg[3]);

		double cut_lj_one = cut_lj_global;
		if (narg == 5) cut_lj_one = force->numeric(FLERR,arg[4]);

		int count = 0;
		for (int i = ilo; i <= ihi; i++) {
		for (int j = MAX(jlo,i); j <= jhi; j++) {
		epsilon[i][j] = epsilon_one;
		sigma[i][j] = sigma_one;
		cut_lj[i][j] = cut_lj_one;
		setflag[i][j] = 1;
		count++;
		}
		}

		if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
		}

		/* ----------------------------------------------------------------------
		init specific to this pair style
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::init_style()
		{
		if (!atom->q_flag)
		error->all(FLERR,"Pair style lj/cut/coul/wolf requires atom attribute q");

		cut_coulsq = cut_coul * cut_coul;

		// request regular or rRESPA neighbor list

		int irequest;
		irequest = neighbor->request(this,instance_me);
		}

		/* ----------------------------------------------------------------------
		init for one type pair i,j and corresponding j,i
		------------------------------------------------------------------------- */

		double PairLJCutCoulWolf::init_one(int i, int j)
		{
		if (setflag[i][j] == 0) {
		epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
		sigma[i][i],sigma[j][j]);
		sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
		cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
		}

		double cut = MAX(cut_lj[i][j],cut_coul);
		cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];

		lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
		lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
		lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
		lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);

		if (offset_flag && (cut_lj[i][j] > 0.0)) {
		double ratio = sigma[i][j] / cut_lj[i][j];
		offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
		} else offset[i][j] = 0.0;

		lj1[j][i] = lj1[i][j];
		lj2[j][i] = lj2[i][j];
		lj3[j][i] = lj3[i][j];
		lj4[j][i] = lj4[i][j];
		offset[j][i] = offset[i][j];

		// compute I,J contribution to long-range tail correction
		// count total # of atoms of type I and J via Allreduce

		if (tail_flag) {
		int *type = atom->type;
		int nlocal = atom->nlocal;

		double count[2],all[2];
		count[0] = count[1] = 0.0;
		for (int k = 0; k < nlocal; k++) {
		if (type[k] == i) count[0] += 1.0;
		if (type[k] == j) count[1] += 1.0;
		}
		MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);

		double sig2 = sigma[i][j]*sigma[i][j];
		double sig6 = sig2sig2sig2;
		double rc3 = cut_lj[i][j]cut_lj[i][j]cut_lj[i][j];
		double rc6 = rc3*rc3;
		double rc9 = rc3*rc6;
		etail_ij = 8.0MY_PIall[0]all[1]epsilon[i][j] *
		sig6 * (sig6 - 3.0rc6) / (9.0rc9);
		ptail_ij = 16.0MY_PIall[0]all[1]epsilon[i][j] *
		sig6 * (2.0sig6 - 3.0rc6) / (9.0*rc9);
		}

		return cut;
		}

		/* ----------------------------------------------------------------------
		proc 0 writes to restart file
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::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(&epsilon[i][j],sizeof(double),1,fp);
		fwrite(&sigma[i][j],sizeof(double),1,fp);
		fwrite(&cut_lj[i][j],sizeof(double),1,fp);
		}
		}
		}

		/* ----------------------------------------------------------------------
		proc 0 reads from restart file, bcasts
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::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(&epsilon[i][j],sizeof(double),1,fp);
		fread(&sigma[i][j],sizeof(double),1,fp);
		fread(&cut_lj[i][j],sizeof(double),1,fp);
		}
		MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
		MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
		MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
		}
		}
		}

		/* ----------------------------------------------------------------------
		proc 0 writes to restart file
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::write_restart_settings(FILE *fp)
		{
		fwrite(&alf,sizeof(double),1,fp);
		fwrite(&cut_lj_global,sizeof(double),1,fp);
		fwrite(&cut_coul,sizeof(double),1,fp);
		fwrite(&offset_flag,sizeof(int),1,fp);
		fwrite(&mix_flag,sizeof(int),1,fp);
		fwrite(&tail_flag,sizeof(int),1,fp);
		}

		/* ----------------------------------------------------------------------
		proc 0 reads from restart file, bcasts
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::read_restart_settings(FILE *fp)
		{
		int me = comm->me;
		if (me == 0) {
		fread(&alf,sizeof(double),1,fp);
		fread(&cut_lj_global,sizeof(double),1,fp);
		fread(&cut_coul,sizeof(double),1,fp);
		fread(&offset_flag,sizeof(int),1,fp);
		fread(&mix_flag,sizeof(int),1,fp);
		fread(&tail_flag,sizeof(int),1,fp);
		}
		MPI_Bcast(&alf,1,MPI_DOUBLE,0,world);
		MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
		MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
		MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
		MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
		MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
		}

		/* ----------------------------------------------------------------------
		proc 0 writes to data file
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::write_data(FILE *fp)
		{
		for (int i = 1; i <= atom->ntypes; i++)
		fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]);
		}

		/* ----------------------------------------------------------------------
		proc 0 writes all pairs to data file
		------------------------------------------------------------------------- */

		void PairLJCutCoulWolf::write_data_all(FILE *fp)
		{
		for (int i = 1; i <= atom->ntypes; i++)
		for (int j = i; j <= atom->ntypes; j++)
		fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]);
		}

src/pair_lj_cut_coul_wolf.h

0 → 100644

+70 −0

File added.

Preview size limit exceeded, changes collapsed.

Admin message