more testing

This commit is contained in:
Steve Plimpton
2022-05-12 16:30:25 -06:00
parent efbfc31868
commit c326bd1caf
9 changed files with 387 additions and 9 deletions

View File

@ -0,0 +1,94 @@
LAMMPS data file for SiO2 zeolite with one methane moleclue
77 atoms
4 atom types
-5.9266 5.9926 xlo xhi
-5.9266 5.9926 ylo yhi
-5.9266 5.9926 zlo zhi
Masses
1 28.0855
2 15.999
3 12.011
4 1.008
Atoms
1 1 1.910418 0.00000 4.38651 2.18123
2 1 1.910418 0.00000 -4.38651 2.18123
3 1 1.910418 0.00000 4.38651 -2.18123
4 1 1.910418 0.00000 -4.38651 -2.18123
5 1 1.910418 2.18123 0.00000 4.38651
6 1 1.910418 2.18123 0.00000 -4.38651
7 1 1.910418 -2.18123 0.00000 4.38651
8 1 1.910418 -2.18123 0.00000 -4.38651
9 1 1.910418 4.38651 2.18123 0.00000
10 1 1.910418 -4.38651 2.18123 0.00000
11 1 1.910418 4.38651 -2.18123 0.00000
12 1 1.910418 -4.38651 -2.18123 0.00000
13 1 1.910418 4.38651 0.00000 -2.18123
14 1 1.910418 -4.38651 0.00000 -2.18123
15 1 1.910418 4.38651 0.00000 2.18123
16 1 1.910418 -4.38651 0.00000 2.18123
17 1 1.910418 0.00000 2.18123 -4.38651
18 1 1.910418 0.00000 2.18123 4.38651
19 1 1.910418 0.00000 -2.18123 -4.38651
20 1 1.910418 0.00000 -2.18123 4.38651
21 1 1.910418 2.18123 4.38651 0.00000
22 1 1.910418 2.18123 -4.38651 0.00000
23 1 1.910418 -2.18123 4.38651 0.00000
24 1 1.910418 -2.18123 -4.38651 0.00000
25 2 0.955209 0.00000 -5.92660 2.64860
26 2 0.955209 0.00000 -5.92660 -2.64860
27 2 0.955209 2.64860 0.00000 -5.92660
28 2 0.955209 -2.64860 0.00000 -5.92660
29 2 0.955209 -5.92660 2.64860 0.00000
30 2 0.955209 -5.92660 -2.64860 0.00000
31 2 0.955209 -5.92660 0.00000 -2.64860
32 2 0.955209 -5.92660 0.00000 2.64860
33 2 0.955209 0.00000 2.64860 -5.92660
34 2 0.955209 0.00000 -2.64860 -5.92660
35 2 0.955209 2.64860 -5.92660 0.00000
36 2 0.955209 -2.64860 -5.92660 0.00000
37 2 0.955209 0.00000 3.45272 3.45272
38 2 0.955209 0.00000 -3.45272 3.45272
39 2 0.955209 0.00000 3.45272 -3.45272
40 2 0.955209 0.00000 -3.45272 -3.45272
41 2 0.955209 3.45272 0.00000 3.45272
42 2 0.955209 3.45272 0.00000 -3.45272
43 2 0.955209 -3.45272 0.00000 3.45272
44 2 0.955209 -3.45272 0.00000 -3.45272
45 2 0.955209 3.45272 3.45272 0.00000
46 2 0.955209 -3.45272 3.45272 0.00000
47 2 0.955209 3.45272 -3.45272 0.00000
48 2 0.955209 -3.45272 -3.45272 0.00000
49 2 0.955209 1.28702 1.28702 4.12598
50 2 0.955209 -1.28702 -1.28702 4.12598
51 2 0.955209 -1.28702 1.28702 -4.12598
52 2 0.955209 1.28702 -1.28702 -4.12598
53 2 0.955209 4.12598 1.28702 1.28702
54 2 0.955209 4.12598 -1.28702 -1.28702
55 2 0.955209 -4.12598 -1.28702 1.28702
56 2 0.955209 -4.12598 1.28702 -1.28702
57 2 0.955209 1.28702 4.12598 1.28702
58 2 0.955209 -1.28702 4.12598 -1.28702
59 2 0.955209 1.28702 -4.12598 -1.28702
60 2 0.955209 -1.28702 -4.12598 1.28702
61 2 0.955209 1.28702 1.28702 -4.12598
62 2 0.955209 -1.28702 -1.28702 -4.12598
63 2 0.955209 1.28702 -1.28702 4.12598
64 2 0.955209 -1.28702 1.28702 4.12598
65 2 0.955209 1.28702 4.12598 -1.28702
66 2 0.955209 -1.28702 4.12598 1.28702
67 2 0.955209 -1.28702 -4.12598 -1.28702
68 2 0.955209 1.28702 -4.12598 1.28702
69 2 0.955209 4.12598 1.28702 -1.28702
70 2 0.955209 4.12598 -1.28702 1.28702
71 2 0.955209 -4.12598 1.28702 1.28702
72 2 0.955209 -4.12598 -1.28702 -1.28702
73 3 -0.66 0.00000 0.00000 0.00000
74 4 0.165 0.00000 -0.89000 -0.62930
75 4 0.165 0.00000 0.89000 -0.62930
76 4 0.165 -0.89000 0.00000 0.62930
77 4 0.165 0.89000 0.00000 0.62930

View File

@ -0,0 +1,46 @@
# MM for water dimer
units real
atom_style full
bond_style harmonic
angle_style harmonic
read_data data.water.dimer
group mm molecule 1
group qm molecule 2
# pair style must define stand-alone short-range Coulombics
pair_style lj/cut/coul/cut 6.0
pair_coeff 1 1 0.13506 3.166
pair_coeff 2 2 0.0 1.0
velocity all create 300.0 458732
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
fix 1 all nve
compute 1 all pair/local dist
compute 2 all reduce max c_1
variable fxabs atom abs(fx)
variable fyabs atom abs(fy)
variable fzabs atom abs(fz)
variable qabs atom abs(q)
compute 3 all reduce max v_fxabs v_fyabs v_fzabs v_qabs
dump 1 all custom 1 dump.water.dimer.mm id x y z q fx fy fz
dump_modify 1 sort id format float "%20.16g"
timestep 0.01
thermo_style custom step cpu temp ke evdwl ecoul epair emol elong pe etotal press &
c_2 c_3[*]
thermo 1
run 1000

View File

@ -0,0 +1,55 @@
# QMMM with NWChem for water dimer
units real
atom_style full
bond_style harmonic
angle_style harmonic
read_data data.water.dimer
group mm molecule 1
group qm molecule 2
# remove bonds/angles in QM water molecule
delete_bonds qm multi remove special
# pair style must define stand-alone short-range Coulombics
pair_style hybrid/overlay lj/cut 6.0 coul/cut 6.0
pair_coeff 1 1 lj/cut 0.13506 3.166
pair_coeff 2 2 lj/cut 0.0 1.0
pair_coeff * * coul/cut
velocity all create 300.0 458732
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
fix 1 all nve
fix 2 qm nwchem template-water.nw water-dimer.nw &
log.pwdft.water-dimer O H
fix_modify 2 energy yes
compute 1 all pair/local dist
compute 2 all reduce max c_1
variable fxabs atom abs(fx)
variable fyabs atom abs(fy)
variable fzabs atom abs(fz)
variable qabs atom abs(q)
compute 3 all reduce max v_fxabs v_fyabs v_fzabs v_qabs
dump 1 all custom 1 dump.water.dimer.qmmm id x y z q fx fy fz
dump_modify 1 sort id format float "%20.16g"
timestep 0.01
thermo_style custom step cpu temp ke evdwl ecoul epair emol elong f_2 pe etotal press &
c_2 c_3[*]
thermo 1
run 1000

View File

@ -0,0 +1,63 @@
# MM for water dimer
units real
atom_style full
bond_style harmonic
angle_style harmonic
read_data data.water.dimer
group mm molecule 1
group qm molecule 2
# pair style must define stand-alone short-range Coulombics
pair_style buck/coul/long 8.0
pair_coeff 1 1 0.13506 3.166
pair_coeff 2 2 0.0 1.0
Si sigma = 0.391995435982 nm
Si epsilon = 2.5104000 kcal/mole ?
qSi = 1.7
qO = -0.85
Si sigma = 3.2
Si eps = 0.0017345 eV
O sigma = 2.7
O eps = 0.02536852 eV
C sigma = 3.4
C eps = 55.055 K
H sigma = 2.65
H eps = 7.901 K
velocity all create 300.0 458732
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes
fix 1 all nve
compute 1 all pair/local dist
compute 2 all reduce max c_1
variable fxabs atom abs(fx)
variable fyabs atom abs(fy)
variable fzabs atom abs(fz)
variable qabs atom abs(q)
compute 3 all reduce max v_fxabs v_fyabs v_fzabs v_qabs
dump 1 all custom 1 dump.water.dimer.mm id x y z q fx fy fz
dump_modify 1 sort id format float "%20.16g"
timestep 0.01
thermo_style custom step cpu temp ke evdwl ecoul epair emol elong pe etotal press &
c_2 c_3[*]
thermo 1
run 1000

42
examples/nwchem/qcoul.py Normal file
View File

@ -0,0 +1,42 @@
from math import sqrt
lines = """
1 1 1 -0.8476 0.13513 -0.05627 0.09445
2 1 2 0.4238 0.79058 0.37197 -0.48187
3 1 2 0.4238 -0.33549 -0.65246 -0.50563
4 2 1 -0.8476 0.02105 0.48666 2.81228
5 2 2 0.4238 -0.73085 1.07786 2.68447
6 2 2 0.4238 0.21349 0.19973 1.90016
"""
qqrd2e = 332.06371
q = []
x = []
lines = lines.split('\n')
for line in lines[1:7]:
print line
id,imol,itype,qone,xone,yone,zone = line.split()
q.append(float(qone))
x.append((float(xone),float(yone),float(zone)))
qp4 = 0.0
qp5 = 0.0
qp6 = 0.0
for i in range(6):
for j in range(6):
if i == j: continue
dx = x[i][0] - x[j][0]
dy = x[i][1] - x[j][1]
dz = x[i][2] - x[j][2]
r = sqrt(dx*dx + dy*dy + dz*dz)
eng = qqrd2e * q[i]*q[j] / r
print "Eng of atoms %d-%d = %g",i+1,j+1,eng
if i+1 == 4: qp4 += eng/q[i]
if i+1 == 5: qp5 += eng/q[i]
if i+1 == 6: qp6 += eng/q[i]
print "QPOT full: 4 %g 5 %g 6 %g" % (qp4,qp5,qp6)

View File

@ -0,0 +1,18 @@
Title "LAMMPS wrapping of PWDFT"
memory 1900 mb
FILEINSERT
echo
GEOMINSERT
nwpw
2d-hcurve
initialize_wavefunction on
cutoff 10.0
xc pbe
end
task pspw gradient

View File

@ -0,0 +1,19 @@
Title "LAMMPS wrapping of PWDFT"
memory 1900 mb
FILEINSERT
echo
GEOMINSERT
nwpw
2d-hcurve
initialize_wavefunction on
cutoff 50.0
mult 1
xc pbe96
end
task pspw gradient

View File

@ -141,6 +141,7 @@ FixNWChem::FixNWChem(LAMMPS *lmp, int narg, char **arg) :
extscalar = 1;
energy_global_flag = 1;
thermo_energy = 1;
comm_forward = 1;
comm_reverse = 1;
// setup unit conversion factors
@ -505,7 +506,6 @@ 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);
@ -525,12 +525,15 @@ void FixNWChem::pre_force_qmmm(int vflag)
}
// reset owned charges to QM values
// communicate changes to ghost atoms
for (int i = 0; i < nqm; i++) {
ilocal = qm2owned[i];
if (ilocal >= 0) q[ilocal] = qqm[i];
}
comm->forward_comm(this);
// reset LAMMPS forces to zero
// NOTE: what about check in force_clear() for external_force_clear = OPENMP ?
// NOTE: what will whichflag be for single snapshot compute of QM forces?
@ -608,9 +611,9 @@ void FixNWChem::post_force_aimd(int vflag)
void FixNWChem::post_force_qmmm(int vflag)
{
int ilocal,jlocal;
double rsq,r2inv,rinv,fpair;
double delta[3];
// int ilocal,jlocal;
// double rsq,r2inv,rinv,fpair;
// double delta[3];
// subtract pairwise QM energy and forces from energy and forces
// LAMMPS just computed for all atoms
@ -621,6 +624,10 @@ void FixNWChem::post_force_qmmm(int vflag)
// this effectively subtract energy from total = pair + kspace + fix
// use xqm & qqm set/computed in pre_force_qmmm() for all QM atoms
// NOTE: no longer needed
// NWChem now subtracts QM/QM contributions from its returned energy and forces
/*
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
@ -675,8 +682,14 @@ void FixNWChem::post_force_qmmm(int vflag)
MPI_Allreduce(&eqm_mine,&eqm,1,MPI_DOUBLE,MPI_SUM,world);
qmenergy -= eqm;
*/
// add NWChem QM forces to owned QM atoms
// do this now, after LAMMPS forces have been re-computed with new QM charges
double **f = atom->f;
int ilocal;
for (int i = 0; i < nqm; i++) {
ilocal = qm2owned[i];
@ -687,11 +700,6 @@ void FixNWChem::post_force_qmmm(int vflag)
}
}
// add QM energy from NWChem
// NOTE: how to calculate this quantity
// trigger per-atom energy computation on next step by pair/kspace
// NOTE: is this needed ?
// only if needed for this fix to calc per-atom forces
@ -730,6 +738,37 @@ void FixNWChem::min_post_force(int vflag)
/* ---------------------------------------------------------------------- */
int FixNWChem::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,m;
double *q = atom->q;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = q[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixNWChem::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
double *q = atom->q;
for (i = first; i < last; i++) q[i] = buf[m++];
}
/* ---------------------------------------------------------------------- */
int FixNWChem::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;

View File

@ -40,6 +40,8 @@ class FixNWChem : public Fix {
void min_post_neighbor() override;
void min_pre_force(int) override;
void min_post_force(int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
double compute_scalar() override;