tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

debugging for QMMM

Browse Source
This commit is contained in:
Steve Plimpton
2022-05-11 15:08:45 -06:00
parent 2410a982c9
commit 7abbb19776
7 changed files with 40 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
examples/nwchem/data.water.dimer
View File

@ -1,14 +1,14 @@
LAMMPS data file for water dimer in large box
6 atoms
4 bonds
2 angles
2 bonds
1 angles
2 atom types
1 bond types
1 angle types
0.0 20.10869277621671 xlo xhi
0.0 20.10869277621671 ylo yhi
0.0 20.10869277621671 zlo zhi
-5.0 5.0 xlo xhi
-5.0 5.0 ylo yhi
-5.0 5.0 zlo zhi
Masses
@ -17,11 +17,11 @@ Masses
Bond Coeffs
1 K r0
1 1.0 1.0
Angle Coeffs
1 K theta0
1 1.0 109.47
Atoms
@ -36,10 +36,7 @@ Bonds
1 1 1 2
2 1 1 3
3 1 4 5
4 1 4 6
Angles
1 1 2 1 3
2 1 5 4 6

2
examples/nwchem/in.w.aimd
View File

@ -15,7 +15,7 @@ neigh_modify delay 0 every 1 check yes
fix 1 all nve
fix 2 all nwchem template-w.nw w.aimd.nw W
fix 2 all nwchem template-w.nw w.aimd.nw log.pwdft.w.aimd W
fix_modify 2 energy yes
dump 1 all custom 1 dump.w.aimd id x y z vx vy vz fx fy fz

3
examples/nwchem/in.w.many
View File

@ -20,7 +20,8 @@ mass 1 183.84
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
fix 2 all nwchem template-w.nw ${datafiles}.nw W
fix 2 all nwchem template-w.nw ${datafiles}.nw &
log.pwdft.${datafiles} W
fix_modify 2 energy yes
dump 1 all custom 1 dump.${datafiles} id x y z fx fy fz

2
examples/nwchem/in.w.single
View File

@ -12,7 +12,7 @@ mass 1 183.84
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
fix 2 all nwchem template-w.nw w.single.nw W
fix 2 all nwchem template-w.nw w.single.nw log.pwdft.w.single W
fix_modify 2 energy yes
dump 1 all custom 1 dump.w.single id x y z fx fy fz

18
examples/nwchem/in.water.dimer
View File

@ -2,19 +2,31 @@
units real
atom_style full
comm_modify cutoff 2.0
bond_style harmonic
angle_style harmonic
read_data data.water.dimer
group mm molecule 1
group qm molecule 2
pair_style hybrid/overlay lj/cut 6.0 coul/cut 6.0
pair_coeff 1 1 lj/cut 1.0 1.0
pair_coeff 2 2 lj/cut 1.0 1.0
pair_coeff * * coul/cut
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
fix 1 all nve
fix 2 all nwchem template.water water-dimer.nw W
fix 2 qm nwchem template-water.nw water-dimer.nw &
SCREEN O H
# log.pwdft.water-dimer O H
fix_modify 2 energy yes
dump 1 all custom 1 dump.water.dimer id x y z vx vy vz fx fy fz
dump 1 all custom 1 dump.water.dimer id x y z fx fy fz
timestep 0.001

3
examples/nwchem/template-water.nw
View File

@ -16,5 +16,4 @@ nwpw
xc pbe96
end
task pspw energy ignore
task pspw optimize ignore
task pspw gradient

24
src/NWCHEM/fix_nwchem.cpp
View File

@ -173,10 +173,10 @@ FixNWChem::FixNWChem(LAMMPS *lmp, int narg, char **arg) :
nqm = ngroup;
if (nqm == atom->natoms) mode = AIMD;
else {
mode = QMMM;
error->all(FLERR,"Fix nwchem for a non-all group (QMMM) is not yet enabled");
}
else mode = QMMM;
//if (mode == QMMM)
// error->all(FLERR,"Fix nwchem for a non-all group (QMMM) is not yet enabled");
memory->create(qmIDs,nqm,"nwchem:qmIDs");
memory->create(xqm,nqm,3,"nwchem:xqm");
@ -347,8 +347,7 @@ void FixNWChem::init()
// one-time initialization of NWChem with its input file
// also one-time setup of qqm = atom charges for all QM atoms
// later calls to NWChem change it for all QM atoms
// changes get copied to LAMMPS per-atom q which is not changed by LAMMPS
// later calls to NWChem change charges for QM atoms
// setup of xqm needed for nwchem_initialize();
if (!qm_init) {
@ -409,7 +408,7 @@ void FixNWChem::pre_force_qmmm(int vflag)
double rsq;
double delta[3];
// invoke pair hybrid sub-style pair coul/long and Kspace directly
// invoke pair hybrid sub-style pair coulomb and Kspace directly
// set eflag = 2 so they calculate per-atom energy
pair_coul->compute(2,0);
@ -450,7 +449,8 @@ void FixNWChem::pre_force_qmmm(int vflag)
// setup 2 NWChem inputs: xqm and qpotential
// xqm = atom coords, mapped into periodic box
// qpotential[i] = Coulomb potential for each atom
// (eatom[i] from pair_coul + kspace) / Qi
// 2 * (eatom[i] from pair_coul + kspace) / Qi
// factor of 2 comes from need to double count energy for each atom
// set for owned atoms, then MPI_Allreduce
// subtract Qj/Rij energy for QM I interacting with all other QM J atoms
// use xqm_mine and qqm_mine for all QM atoms
@ -466,7 +466,7 @@ void FixNWChem::pre_force_qmmm(int vflag)
ilocal = qm2owned[i];
if (ilocal >= 0) {
if (q[ilocal] == 0.0) qpotential_mine[i] = 0.0;
else qpotential_mine[i] = ecoul[ilocal] / q[ilocal];
else qpotential_mine[i] = 2.0 * ecoul[ilocal] / q[ilocal];
}
}
@ -505,8 +505,10 @@ void FixNWChem::pre_force_qmmm(int vflag)
// fqm,qqm = forces & charges
// qmenergy = QM energy of entire system
int nwerr = lammps_pspw_qmmm_minimizer(world,&xqm[0][0],qpotential,
&fqm[0][0],qqm,&qmenergy);
//int nwerr = dummy_pspw_qmmm_minimizer(world,&xqm[0][0],qpotential,
// &fqm[0][0],qqm,&qmenergy);
@ -575,7 +577,6 @@ void FixNWChem::post_force_aimd(int vflag)
// qmenergy = QM energy of entire system
int nwerr = lammps_pspw_aimd_minimizer(world,&xqm[0][0],&fqm[0][0],&qmenergy);
if (nwerr) error->all(FLERR,"Internal NWChem error");
// unit conversion from NWChem to LAMMPS
@ -813,7 +814,8 @@ void FixNWChem::nwchem_initialize()
if (strcmp(line,"FILEINSERT\n") == 0) {
fprintf(infile,"start %s\n",nw_input);
if (strcmp(nw_output,"NULL") == 0) fprintf(infile,"print off\n");
if (strcmp(nw_output,"SCREEN") == 0) {
} else if (strcmp(nw_output,"NULL") == 0) fprintf(infile,"print off\n");
else fprintf(infile,"redirect_filename %s\n",nw_output);
continue;
}