From 43b81297d7b256f1527aeb667ad871816f122dcc Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 6 Jan 2012 18:31:54 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7448 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-EFF/README | 15 ++ src/USER-EFF/compute_ke_atom_eff.cpp | 4 +- src/USER-EFF/compute_ke_eff.cpp | 3 +- src/USER-EFF/compute_temp_deform_eff.cpp | 13 +- src/USER-EFF/compute_temp_eff.cpp | 12 +- src/USER-EFF/compute_temp_region_eff.cpp | 23 +- src/USER-EFF/fix_langevin_eff.cpp | 280 +++++++++++++++++++---- src/USER-EFF/fix_nh_eff.cpp | 4 +- src/USER-EFF/fix_nve_eff.cpp | 7 +- src/USER-EFF/pair_eff_cut.cpp | 12 +- 10 files changed, 298 insertions(+), 75 deletions(-) diff --git a/src/USER-EFF/README b/src/USER-EFF/README index 6bae419b1c..aaa11c3494 100644 --- a/src/USER-EFF/README +++ b/src/USER-EFF/README @@ -74,3 +74,18 @@ real units, upgraded post-processing and visualization tools, added support for "compute pair eff" to allow thermo prints with the different eFF energy components (eke, epauli, ecoul and errestrain), fixed radial scaling factors in the eff langevin thermostat. + +12/2011: Added support for "zero" option in fix langevin/eff (see doc), and +adjusted fix_langevin_eff.cpp to correctly thermostat between nuclear and electronic dof +(required additional scaling of friction term in the Langevin equations of motion). +Radial electron mass now scales as a function of system dimension. + +Bug fixes: +(10-2011): Thanks to Christian Chenard-Lemire (U Montreal) for reporting a bug in the +fixed pair_eff_cut.cpp fixedcore-pseudocore interactions (an incorrect index and a missing +elec-core call to account for the 2 electrons from the fixed core.) +(12-2011): Corrected undefined natom variable in fix_langevin_eff (recent changes in +main fix_langevin class caused compilation error in user-eff). +Corrected thermostat in fix langevin/eff as described in version 12/2011. + + diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp index 1b68b60d22..695cdc8874 100644 --- a/src/USER-EFF/compute_ke_atom_eff.cpp +++ b/src/USER-EFF/compute_ke_atom_eff.cpp @@ -25,6 +25,7 @@ #include "force.h" #include "memory.h" #include "error.h" +#include "domain.h" using namespace LAMMPS_NS; @@ -90,6 +91,7 @@ void ComputeKEAtomEff::compute_peratom() int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; if (mass) for (int i = 0; i < nlocal; i++) { @@ -97,7 +99,7 @@ void ComputeKEAtomEff::compute_peratom() ke[i] = 0.5 * mvv2e * mass[type[i]] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); if (fabs(spin[i])==1) - ke[i] += 0.5 * mvv2e * mass[type[i]] * ervel[i]*ervel[i] * 0.75; + ke[i] += 0.5 * mvv2e * mefactor * mass[type[i]] * ervel[i]*ervel[i]; } else ke[i] = 0.0; } } diff --git a/src/USER-EFF/compute_ke_eff.cpp b/src/USER-EFF/compute_ke_eff.cpp index 43a2a952f8..5da7d25740 100644 --- a/src/USER-EFF/compute_ke_eff.cpp +++ b/src/USER-EFF/compute_ke_eff.cpp @@ -63,6 +63,7 @@ double ComputeKEEff::compute_scalar() int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; double ke = 0.0; @@ -71,7 +72,7 @@ double ComputeKEEff::compute_scalar() if (mask[i] & groupbit) { ke += mass[type[i]]*(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); - if (fabs(spin[i])==1) ke += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) ke += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp index ccba4d782c..68e73a83d1 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -113,6 +113,7 @@ void ComputeTempDeformEff::dof_compute() int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); + // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); @@ -135,6 +136,7 @@ double ComputeTempDeformEff::compute_scalar() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; // lamda = 0-1 triclinic lamda coords // vstream = streaming velocity = Hrate*lamda + Hratelo @@ -159,7 +161,7 @@ double ComputeTempDeformEff::compute_scalar() if (mass) { t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2])* mass[type[i]]; - if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } @@ -185,7 +187,8 @@ void ComputeTempDeformEff::compute_vector() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - + double mefactor = domain->dimension/4.0; + double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; @@ -211,9 +214,9 @@ void ComputeTempDeformEff::compute_vector() t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; if (fabs(spin[i])==1) { - t[0] += 0.75 * massone * ervel[i]*ervel[i]; - t[1] += 0.75 * massone * ervel[i]*ervel[i]; - t[2] += 0.75 * massone * ervel[i]*ervel[i]; + t[0] += mefactor * massone * ervel[i]*ervel[i]; + t[1] += mefactor * massone * ervel[i]*ervel[i]; + t[2] += mefactor * massone * ervel[i]*ervel[i]; } } diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index e2089a0d60..cf406c53e6 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -84,7 +84,7 @@ void ComputeTempEff::dof_compute() int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); - // average over nuclear dof only + // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; @@ -105,6 +105,7 @@ double ComputeTempEff::compute_scalar() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; double t = 0.0; @@ -113,7 +114,7 @@ double ComputeTempEff::compute_scalar() if (mask[i] & groupbit) { t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; - if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } } @@ -139,6 +140,7 @@ void ComputeTempEff::compute_vector() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; @@ -153,9 +155,9 @@ void ComputeTempEff::compute_vector() t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; if (fabs(spin[i])==1) { - t[0] += 0.75*massone*ervel[i]*ervel[i]; - t[1] += 0.75*massone*ervel[i]*ervel[i]; - t[2] += 0.75*massone*ervel[i]*ervel[i]; + t[0] += mefactor*massone*ervel[i]*ervel[i]; + t[1] += mefactor*massone*ervel[i]*ervel[i]; + t[2] += mefactor*massone*ervel[i]*ervel[i]; } } diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp index f43e92095f..0dd2bf3cf6 100644 --- a/src/USER-EFF/compute_temp_region_eff.cpp +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -104,9 +104,11 @@ double ComputeTempRegionEff::compute_scalar() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; Region *region = domain->regions[iregion]; int count = 0; + int ecount = 0; double t = 0.0; if (mass) { @@ -115,12 +117,16 @@ double ComputeTempRegionEff::compute_scalar() count++; t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; - if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) { + t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; + ecount++; + } } } double tarray[2],tarray_all[2]; - tarray[0] = count; + // Assume 3/2 k T per nucleus + tarray[0] = count-ecount; tarray[1] = t; MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world); dof = domain->dimension * tarray_all[0] - extra_dof; @@ -130,12 +136,6 @@ double ComputeTempRegionEff::compute_scalar() if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (fabs(spin[i])==1) one++; } - int nelectrons_region; - MPI_Allreduce(&one,&nelectrons_region,1,MPI_INT,MPI_SUM,world); - - // average over nuclear dof only - - dof -= domain->dimension * nelectrons_region ; if (dof > 0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); else scalar = 0.0; @@ -158,6 +158,7 @@ void ComputeTempRegionEff::compute_vector() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; Region *region = domain->regions[iregion]; double massone,t[6]; @@ -175,9 +176,9 @@ void ComputeTempRegionEff::compute_vector() t[5] += massone * v[i][1]*v[i][2]; if (fabs(spin[i])==1) { - t[0] += 0.75 * massone * ervel[i]*ervel[i]; - t[1] += 0.75 * massone * ervel[i]*ervel[i]; - t[2] += 0.75 * massone * ervel[i]*ervel[i]; + t[0] += mefactor * massone * ervel[i]*ervel[i]; + t[1] += mefactor * massone * ervel[i]*ervel[i]; + t[2] += mefactor * massone * ervel[i]*ervel[i]; } } diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp index 7328fdb23b..faff9909d1 100644 --- a/src/USER-EFF/fix_langevin_eff.cpp +++ b/src/USER-EFF/fix_langevin_eff.cpp @@ -20,6 +20,7 @@ #include "string.h" #include "stdlib.h" #include "fix_langevin_eff.h" +#include "math_extra.h" #include "atom.h" #include "force.h" #include "update.h" @@ -29,13 +30,20 @@ #include "region.h" #include "respa.h" #include "comm.h" +#include "input.h" +#include "variable.h" #include "random_mars.h" #include "memory.h" #include "error.h" +#include "group.h" using namespace LAMMPS_NS; enum{NOBIAS,BIAS}; +enum{CONSTANT,EQUAL,ATOM}; + +#define SINERTIA 0.4 // moment of inertia prefactor for sphere +#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid /* ---------------------------------------------------------------------- */ @@ -56,22 +64,51 @@ FixLangevinEff::~FixLangevinEff() void FixLangevinEff::post_force_no_tally() { - double gamma1,gamma2; + double gamma1,gamma2,t_target; double **v = atom->v; double **f = atom->f; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - double *ervel = atom->ervel; double *erforce = atom->erforce; int *spin = atom->spin; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; + double sqrtmefactor = sqrt(mefactor); double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; - double t_target = t_start + delta * (t_stop-t_start); - double tsqrt = sqrt(t_target); + + // set current t_target and t_sqrt + // if variable temp, evaluate variable, wrap with clear/add + // reallocate tforce array if necessary + + if (tstyle == CONSTANT) { + t_target = t_start + delta * (t_stop-t_start); + tsqrt = sqrt(t_target); + } else { + modify->clearstep_compute(); + if (tstyle == EQUAL) { + t_target = input->variable->compute_equal(tvar); + if (t_target < 0.0) + error->one(FLERR,"Fix langevin/eff variable returned negative temperature"); + tsqrt = sqrt(t_target); + } else { + if (nlocal > maxatom2) { + maxatom2 = atom->nmax; + memory->destroy(tforce); + memory->create(tforce,maxatom2,"langevin/eff:tforce"); + } + 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/eff variable returned negative temperature"); + } + modify->addstep_compute(update->ntimestep + 1); + } // apply damping and thermostat to atoms in group // for BIAS: @@ -79,72 +116,179 @@ void FixLangevinEff::post_force_no_tally() // computed on current nlocal atoms to remove bias // test v = 0 since some computes mask non-participating atoms via v = 0 // and added force has extra term not multiplied by v = 0 + // for ZEROFLAG: + // sum random force over all atoms in group + // subtract sum/particles from each atom in group + + double fran[4],fsum[4],fsumall[4]; + fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; + + double boltz = force->boltz; + double dt = update->dt; + double mvv2e = force->mvv2e; + double ftm2v = force->ftm2v; + + int particles = group->count(igroup); + if (zeroflag) { + if (particles == 0) + error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons"); + } + + // find number of electrons in group + int dof,fix_dof; + dof = domain->dimension * particles; + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + + // extra_dof = domain->dimension + dof -= domain->dimension + fix_dof; + + int one = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (fabs(spin[i])==1) one++; + } + int nelectrons, dofelectrons, dofnuclei; + MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); + dofelectrons = domain->dimension*nelectrons; + dofnuclei = dof-dofelectrons; + + // thermal partitioning factor between nuclei and electrons + // extra dof from electron size + double gfactor3=(double) (dof+nelectrons)/dofnuclei; if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; gamma2 = gfactor2[type[i]] * tsqrt; - f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5); - f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5); - f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - if (abs(spin[i])==1) - erforce[i] += 0.75*gamma1*ervel[i] + - 0.866025404*gamma2*(random->uniform()-0.5); + fran[0] = gamma2*(random->uniform()-0.5); + fran[1] = gamma2*(random->uniform()-0.5); + fran[2] = gamma2*(random->uniform()-0.5); + f[i][0] += gamma1*v[i][0] + fran[0]; + f[i][1] += gamma1*v[i][1] + fran[1]; + f[i][2] += gamma1*v[i][2] + fran[2]; + fsum[0] += fran[0]; + fsum[1] += fran[1]; + fsum[2] += fran[2]; + if (fabs(spin[i])==1) { + fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5); + erforce[i] += mefactor*gamma1*ervel[i]+fran[3]; + fsum[3] += fran[3]; + } } } } else if (which == BIAS) { - double tmp = temperature->compute_scalar(); + temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; gamma2 = gfactor2[type[i]] * tsqrt; temperature->remove_bias(i,v[i]); + fran[0] = gamma2*(random->uniform()-0.5); + fran[1] = gamma2*(random->uniform()-0.5); + fran[2] = gamma2*(random->uniform()-0.5); if (v[i][0] != 0.0) - f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5); + f[i][0] += gamma1*v[i][0] + fran[0]; if (v[i][1] != 0.0) - f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5); + f[i][1] += gamma1*v[i][1] + fran[1]; if (v[i][2] != 0.0) - f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - if (abs(spin[i])==1 && ervel[i] != 0.0) - erforce[i] += 0.75*gamma1*ervel[i] + - 0.866025404*gamma2*(random->uniform()-0.5); + f[i][2] += gamma1*v[i][2] + fran[2]; + fsum[0] += fran[0]; + fsum[1] += fran[1]; + fsum[2] += fran[2]; + if (fabs(spin[i])==1) { + fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5); + if (ervel[i] != 0.0) erforce[i] += mefactor*gamma1*ervel[i]+fran[3]; + fsum[3] += fran[3]; + } temperature->restore_bias(i,v[i]); } } } + + // set total force to zero + + if (zeroflag) { + MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world); + fsumall[0] /= particles; + fsumall[1] /= particles; + fsumall[2] /= particles; + fsumall[3] /= nelectrons; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + f[i][0] -= fsumall[0]; + f[i][1] -= fsumall[1]; + f[i][2] -= fsumall[2]; + if (fabs(spin[i])==1) erforce[i] -= fsumall[3]; + } + } + } } /* ---------------------------------------------------------------------- */ void FixLangevinEff::post_force_tally() { - double gamma1,gamma2; + double gamma1,gamma2,t_target; - // reallocate flangevin if necessary + // reallocate flangevin and erforcelangevin if necessary - if (atom->nmax > nmax) { + if (atom->nlocal > maxatom1) { memory->destroy(flangevin); memory->destroy(erforcelangevin); - nmax = atom->nmax; - memory->create(flangevin,nmax,3,"langevin:flangevin"); - memory->create(erforcelangevin,nmax,"langevin/eff:erforcelangevin"); + maxatom1 = atom->nmax; + memory->create(flangevin,maxatom1,3,"langevin/eff:flangevin"); + memory->create(erforcelangevin,maxatom1,"langevin/eff:erforcelangevin"); } double **v = atom->v; double **f = atom->f; + double *erforce = atom->erforce; + double *ervel = atom->ervel; + int *spin = atom->spin; + double mefactor = domain->dimension/4.0; + double sqrtmefactor = sqrt(mefactor); + int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - double *erforce = atom->erforce; - double *ervel = atom->ervel; - int *spin = atom->spin; - double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; - double t_target = t_start + delta * (t_stop-t_start); - double tsqrt = sqrt(t_target); + + // set current t_target and t_sqrt + // if variable temp, evaluate variable, wrap with clear/add + // reallocate tforce array if necessary + + if (tstyle == CONSTANT) { + t_target = t_start + delta * (t_stop-t_start); + tsqrt = sqrt(t_target); + } else { + modify->clearstep_compute(); + if (tstyle == EQUAL) { + t_target = input->variable->compute_equal(tvar); + if (t_target < 0.0) + error->one(FLERR,"Fix langevin/eff variable returned negative temperature"); + tsqrt = sqrt(t_target); + } else { + if (nlocal > maxatom2) { + maxatom2 = atom->nmax; + memory->destroy(tforce); + memory->create(tforce,maxatom2,"langevin/eff:tforce"); + } + 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/eff variable returned negative temperature"); + } + modify->addstep_compute(update->ntimestep + 1); + } // apply damping and thermostat to appropriate atoms // for BIAS: @@ -153,42 +297,81 @@ void FixLangevinEff::post_force_tally() // test v = 0 since some computes mask non-participating atoms via v = 0 // and added force has extra term not multiplied by v = 0 + double boltz = force->boltz; + double dt = update->dt; + double mvv2e = force->mvv2e; + double ftm2v = force->ftm2v; + + int particles = group->count(igroup); + if (zeroflag) { + if (particles == 0) + error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons"); + } + + // find number of electrons in group + int dof,fix_dof; + dof = domain->dimension * particles; + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + + // extra_dof = domain->dimension + dof -= domain->dimension + fix_dof; + + int one = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (fabs(spin[i])==1) one++; + } + int nelectrons, dofelectrons, dofnuclei; + MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); + dofelectrons = domain->dimension*nelectrons; + dofnuclei = dof-dofelectrons; + + // thermal partitioning factor between nuclei and electrons + // extra dof from electron size + double gfactor3=(double) (dof+nelectrons)/dofnuclei; + if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; gamma2 = gfactor2[type[i]] * tsqrt; flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5); flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5); flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - erforcelangevin[i] = 0.75*gamma1*ervel[i] + - 0.866025404*gamma2*(random->uniform()-0.5); f[i][0] += flangevin[i][0]; f[i][1] += flangevin[i][1]; f[i][2] += flangevin[i][2]; - if (abs(spin[i])==1) erforce[i] += erforcelangevin[i]; + if (fabs(spin[i])==1) { + erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5); + erforce[i] += erforcelangevin[i]; + } } } } else if (which == BIAS) { - double tmp = temperature->compute_scalar(); + temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; gamma2 = gfactor2[type[i]] * tsqrt; temperature->remove_bias(i,v[i]); flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5); flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5); flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - erforcelangevin[i] = 0.75*gamma1*ervel[i] + - 0.866025404*gamma2*(random->uniform()-0.5); if (v[i][0] != 0.0) f[i][0] += flangevin[i][0]; else flangevin[i][0] = 0.0; if (v[i][1] != 0.0) f[i][1] += flangevin[i][1]; else flangevin[i][1] = 0.0; if (v[i][2] != 0.0) f[i][2] += flangevin[i][2]; else flangevin[i][2] = 0.0; - if (abs(spin[i])==1 && ervel[i] != 0.0) - erforce[i] += erforcelangevin[i]; + if (fabs(spin[i])==1) { + erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5); + if (ervel[i] != 0.0) erforce[i] += erforcelangevin[i]; + else erforcelangevin[i] = 0.0; + } temperature->restore_bias(i,v[i]); } } @@ -213,8 +396,8 @@ void FixLangevinEff::end_of_step() 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 (abs(spin[i])==1) energy_onestep += erforcelangevin[i]; + flangevin[i][2]*v[i][2]; + if (fabs(spin[i])==1) energy_onestep += erforcelangevin[i]; } energy += energy_onestep*update->dt; } @@ -223,7 +406,7 @@ void FixLangevinEff::end_of_step() double FixLangevinEff::compute_scalar() { - if (!tally) return 0.0; + if (!tally || flangevin == NULL || erforcelangevin == NULL) return 0.0; // capture the very first energy transfer to thermal reservoir @@ -238,7 +421,7 @@ double FixLangevinEff::compute_scalar() 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 (abs(spin[i])==1) energy_onestep += erforcelangevin[i]; + if (fabs(spin[i])==1) energy_onestep += erforcelangevin[i]; } energy = 0.5*energy_onestep*update->dt; } @@ -249,3 +432,4 @@ double FixLangevinEff::compute_scalar() MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); return -energy_all; } + diff --git a/src/USER-EFF/fix_nh_eff.cpp b/src/USER-EFF/fix_nh_eff.cpp index e3dde84178..12bcb90760 100644 --- a/src/USER-EFF/fix_nh_eff.cpp +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -21,6 +21,7 @@ #include "atom_vec.h" #include "group.h" #include "error.h" +#include "domain.h" using namespace LAMMPS_NS; @@ -48,6 +49,7 @@ void FixNHEff::nve_v() double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; + double mefactor = domain->dimension/4.0; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -60,7 +62,7 @@ void FixNHEff::nve_v() if (mask[i] & groupbit) { if (fabs(spin[i])==1) { dtfm = dtf / mass[type[i]]; - ervel[i] = dtfm * erforce[i] / 0.75; + ervel[i] = dtfm * erforce[i] / mefactor; } } } diff --git a/src/USER-EFF/fix_nve_eff.cpp b/src/USER-EFF/fix_nve_eff.cpp index 1466c777df..9b4613df3e 100644 --- a/src/USER-EFF/fix_nve_eff.cpp +++ b/src/USER-EFF/fix_nve_eff.cpp @@ -25,6 +25,7 @@ #include "respa.h" #include "error.h" #include "math.h" +#include "domain.h" using namespace LAMMPS_NS; @@ -80,6 +81,7 @@ void FixNVEEff::initial_integrate(int vflag) double *erforce = atom->erforce; double *mass = atom->mass; int *spin = atom->spin; + double mefactor = domain->dimension/4.0; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -98,7 +100,7 @@ void FixNVEEff::initial_integrate(int vflag) x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; if (fabs(spin[i])==1) { - ervel[i] += dtfm * erforce[i] / 0.75; + ervel[i] += dtfm * erforce[i] / mefactor; eradius[i] += dtv * ervel[i]; } } @@ -118,6 +120,7 @@ void FixNVEEff::final_integrate() double **f = atom->f; double *mass = atom->mass; int *spin = atom->spin; + double mefactor = domain->dimension/4.0; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -133,7 +136,7 @@ void FixNVEEff::final_integrate() v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; if (fabs(spin[i])==1) - ervel[i] += dtfm * erforce[i] / 0.75; + ervel[i] += dtfm * erforce[i] / mefactor; } } } diff --git a/src/USER-EFF/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp index 8f99134239..ee191a0d1d 100644 --- a/src/USER-EFF/pair_eff_cut.cpp +++ b/src/USER-EFF/pair_eff_cut.cpp @@ -402,6 +402,11 @@ void PairEffCut::compute(int eflag, int vflag) ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair); ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul, &fpair,&e2rforce); + ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul, + &fpair,&e2rforce); + PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair, + &ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B, + PAULI_CORE_C); PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair, &ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B, PAULI_CORE_C); @@ -427,7 +432,7 @@ void PairEffCut::compute(int eflag, int vflag) else if ((abs(spin[i]) == 1 || spin[i] == 2) && spin[j] == 3) { e1rforce = ecp_e1rforce = 0.0; - if (abs(spin[j]) == 1) { + if (abs(spin[i]) == 1) { ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul, &fpair,&e1rforce); PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair, @@ -438,6 +443,11 @@ void PairEffCut::compute(int eflag, int vflag) ElecCoreNuc(qxq,rc,eradius[i],&ecoul,&fpair); ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul, &fpair,&e1rforce); + ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul, + &fpair,&e1rforce); + PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair, + &ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B, + PAULI_CORE_C); PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair, &ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B, PAULI_CORE_C);