diff --git a/examples/QM/LATTE/2uo2.lmp b/examples/QUANTUM/LATTE/2uo2.lmp similarity index 100% rename from examples/QM/LATTE/2uo2.lmp rename to examples/QUANTUM/LATTE/2uo2.lmp diff --git a/examples/QM/LATTE/3uo2.lmp b/examples/QUANTUM/LATTE/3uo2.lmp similarity index 100% rename from examples/QM/LATTE/3uo2.lmp rename to examples/QUANTUM/LATTE/3uo2.lmp diff --git a/examples/QM/LATTE/4uo2.lmp b/examples/QUANTUM/LATTE/4uo2.lmp similarity index 100% rename from examples/QM/LATTE/4uo2.lmp rename to examples/QUANTUM/LATTE/4uo2.lmp diff --git a/examples/QM/LATTE/README b/examples/QUANTUM/LATTE/README similarity index 100% rename from examples/QM/LATTE/README rename to examples/QUANTUM/LATTE/README diff --git a/examples/QM/LATTE/bondints.table b/examples/QUANTUM/LATTE/bondints.table similarity index 100% rename from examples/QM/LATTE/bondints.table rename to examples/QUANTUM/LATTE/bondints.table diff --git a/examples/QM/LATTE/dump.8Sep22.aimd.mpi.1 b/examples/QUANTUM/LATTE/dump.8Sep22.aimd.mpi.1 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.aimd.mpi.1 rename to examples/QUANTUM/LATTE/dump.8Sep22.aimd.mpi.1 diff --git a/examples/QM/LATTE/dump.8Sep22.aimd.mpi.2 b/examples/QUANTUM/LATTE/dump.8Sep22.aimd.mpi.2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.aimd.mpi.2 rename to examples/QUANTUM/LATTE/dump.8Sep22.aimd.mpi.2 diff --git a/examples/QM/LATTE/dump.8Sep22.aimd.plugin b/examples/QUANTUM/LATTE/dump.8Sep22.aimd.plugin similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.aimd.plugin rename to examples/QUANTUM/LATTE/dump.8Sep22.aimd.plugin diff --git a/examples/QM/LATTE/dump.8Sep22.series.mpi.1.2uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.1.2uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.mpi.1.2uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.1.2uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.mpi.1.3uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.1.3uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.mpi.1.3uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.1.3uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.mpi.1.4uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.1.4uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.mpi.1.4uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.1.4uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.mpi.2.2uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.2.2uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.mpi.2.2uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.2.2uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.mpi.2.3uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.2.3uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.mpi.2.3uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.2.3uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.mpi.2.4uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.2.4uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.mpi.2.4uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.mpi.2.4uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.plugin.2uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.plugin.2uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.plugin.2uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.plugin.2uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.plugin.3uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.plugin.3uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.plugin.3uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.plugin.3uo2 diff --git a/examples/QM/LATTE/dump.8Sep22.series.plugin.4uo2 b/examples/QUANTUM/LATTE/dump.8Sep22.series.plugin.4uo2 similarity index 100% rename from examples/QM/LATTE/dump.8Sep22.series.plugin.4uo2 rename to examples/QUANTUM/LATTE/dump.8Sep22.series.plugin.4uo2 diff --git a/examples/QM/LATTE/electrons.dat b/examples/QUANTUM/LATTE/electrons.dat similarity index 100% rename from examples/QM/LATTE/electrons.dat rename to examples/QUANTUM/LATTE/electrons.dat diff --git a/examples/QM/LATTE/in.aimd b/examples/QUANTUM/LATTE/in.aimd similarity index 100% rename from examples/QM/LATTE/in.aimd rename to examples/QUANTUM/LATTE/in.aimd diff --git a/examples/QM/LATTE/in.aimd.plugin b/examples/QUANTUM/LATTE/in.aimd.plugin similarity index 100% rename from examples/QM/LATTE/in.aimd.plugin rename to examples/QUANTUM/LATTE/in.aimd.plugin diff --git a/examples/QM/LATTE/in.series b/examples/QUANTUM/LATTE/in.series similarity index 100% rename from examples/QM/LATTE/in.series rename to examples/QUANTUM/LATTE/in.series diff --git a/examples/QM/LATTE/in.series.plugin b/examples/QUANTUM/LATTE/in.series.plugin similarity index 100% rename from examples/QM/LATTE/in.series.plugin rename to examples/QUANTUM/LATTE/in.series.plugin diff --git a/examples/QM/LATTE/latte.in b/examples/QUANTUM/LATTE/latte.in similarity index 100% rename from examples/QM/LATTE/latte.in rename to examples/QUANTUM/LATTE/latte.in diff --git a/examples/QM/LATTE/log.8Sep22.aimd.lammps.mpi.1 b/examples/QUANTUM/LATTE/log.8Sep22.aimd.lammps.mpi.1 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.aimd.lammps.mpi.1 rename to examples/QUANTUM/LATTE/log.8Sep22.aimd.lammps.mpi.1 diff --git a/examples/QM/LATTE/log.8Sep22.aimd.lammps.mpi.2 b/examples/QUANTUM/LATTE/log.8Sep22.aimd.lammps.mpi.2 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.aimd.lammps.mpi.2 rename to examples/QUANTUM/LATTE/log.8Sep22.aimd.lammps.mpi.2 diff --git a/examples/QM/LATTE/log.8Sep22.aimd.lammps.plugin b/examples/QUANTUM/LATTE/log.8Sep22.aimd.lammps.plugin similarity index 100% rename from examples/QM/LATTE/log.8Sep22.aimd.lammps.plugin rename to examples/QUANTUM/LATTE/log.8Sep22.aimd.lammps.plugin diff --git a/examples/QM/LATTE/log.8Sep22.aimd.latte.1 b/examples/QUANTUM/LATTE/log.8Sep22.aimd.latte.1 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.aimd.latte.1 rename to examples/QUANTUM/LATTE/log.8Sep22.aimd.latte.1 diff --git a/examples/QM/LATTE/log.8Sep22.aimd.latte.2 b/examples/QUANTUM/LATTE/log.8Sep22.aimd.latte.2 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.aimd.latte.2 rename to examples/QUANTUM/LATTE/log.8Sep22.aimd.latte.2 diff --git a/examples/QM/LATTE/log.8Sep22.aimd.latte.plugin b/examples/QUANTUM/LATTE/log.8Sep22.aimd.latte.plugin similarity index 100% rename from examples/QM/LATTE/log.8Sep22.aimd.latte.plugin rename to examples/QUANTUM/LATTE/log.8Sep22.aimd.latte.plugin diff --git a/examples/QM/LATTE/log.8Sep22.series.lammps.mpi.1 b/examples/QUANTUM/LATTE/log.8Sep22.series.lammps.mpi.1 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.series.lammps.mpi.1 rename to examples/QUANTUM/LATTE/log.8Sep22.series.lammps.mpi.1 diff --git a/examples/QM/LATTE/log.8Sep22.series.lammps.mpi.2 b/examples/QUANTUM/LATTE/log.8Sep22.series.lammps.mpi.2 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.series.lammps.mpi.2 rename to examples/QUANTUM/LATTE/log.8Sep22.series.lammps.mpi.2 diff --git a/examples/QM/LATTE/log.8Sep22.series.lammps.plugin b/examples/QUANTUM/LATTE/log.8Sep22.series.lammps.plugin similarity index 100% rename from examples/QM/LATTE/log.8Sep22.series.lammps.plugin rename to examples/QUANTUM/LATTE/log.8Sep22.series.lammps.plugin diff --git a/examples/QM/LATTE/log.8Sep22.series.latte.1 b/examples/QUANTUM/LATTE/log.8Sep22.series.latte.1 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.series.latte.1 rename to examples/QUANTUM/LATTE/log.8Sep22.series.latte.1 diff --git a/examples/QM/LATTE/log.8Sep22.series.latte.2 b/examples/QUANTUM/LATTE/log.8Sep22.series.latte.2 similarity index 100% rename from examples/QM/LATTE/log.8Sep22.series.latte.2 rename to examples/QUANTUM/LATTE/log.8Sep22.series.latte.2 diff --git a/examples/QM/LATTE/ppots.dftb b/examples/QUANTUM/LATTE/ppots.dftb similarity index 100% rename from examples/QM/LATTE/ppots.dftb rename to examples/QUANTUM/LATTE/ppots.dftb diff --git a/examples/QUANTUM/PySCF/README b/examples/QUANTUM/PySCF/README new file mode 100644 index 0000000000..775ceb8015 --- /dev/null +++ b/examples/QUANTUM/PySCF/README @@ -0,0 +1,103 @@ +# Test runs of QMMM with LAMMPS and PySCF + +Step 1: build LAMMPS +Step 2: download/build the MDI code coupling package +Step 3: download/build or install PySCF +Step 4: run 2-water QMMM problem for a few steps + +--------------------------------- +--------------------------------- + +Step 1: build LAMMPS + +Traditional make: + +% cd ~/lammps/src +% make yes-mdi yes-molecule yes-kspace +% make -j mpi +% cp lmp_mpi ~/lammps/examples/QUANTUM/PySCF + +CMake: + +% cd ~/lammps +% mkdir build_test; cd build_test +% cmake ../cmake -D PKG_MDI=yes -D PKG_MOLECULE=yes -D PKG_KSPACE=yes +% make -j +% cp lmp_mpi ~/lammps/examples/QUANTUM/PySCF + +Note: the molecule and kspace packages are needed for the water test +problem. + +--------------------------------- +--------------------------------- + +Step 2: download/build the MDI code coupling package + +(a) grab the MDI Git repo + +% mkdir git; cd mdi +% git clone git@github.com:MolSSI-MDI/MDI_Library.git git + +(b) build MDI + +% cd mdi/git +% mkdir build; cd build +% cmake .. # include support for all langauges +% make -j + +(c) install mdi.py into your Python: + +% cd mdi/git +% pip3 install . + +--------------------------------- +--------------------------------- + +Step 3: download/build or install PySCF + +(a) install PySCF on your box + +NOTE: add instructions here + +(b) Add something similar to the following to your .bashrc or .cshrc +file so that Python can find PySCF: + +For bash: + +% export PYTHONPATH="$PYTHONPATH:/home/sjplimp/pyscf/git" +% hash -r + +For (t)csh: + +% setenv PYTHONPATH ${PYTHONPATH}:/home/sjplimp/pyscf/git +% rehash + +(c) Check that you can import the 4 Python modules which the script +that wraps PySCF will need: + +% python +>>> import numpy as np +>>> from mpi4py import MPI +>>> import MDI_Library as mdi +>>> import pyscf + +--------------------------------- +--------------------------------- + +Step 4: run 2-water QMMM problem for a few steps + +% cd ~/lammps/examples/QUANTUM/PySCF + +# Run with TCP: 1 proc each + +% lmp_mpi -mdi "-name LMP -role DRIVER -method TCP -port 8021" -log log.water.pyscf.qmmm.tcp.1 -in in.water.pyscf.qmmm & + +% python pyscf_mdi.py -mdi "-name PYSCF -role ENGINE -method TCP -port 8021 -hostname localhost" + +# Run with MPI: 1 proc each + +% mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.water.pyscf.qmmm.mpi.1 -in in.water.pyscf.qmmm : -np 1 python pyscf_mdi.py -mdi "-name PYSCF -role ENGINE -method MPI" -pbc yes + +# Run in plugin mode: 1 proc + +% lmp_mpi -mdi "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/work_qm" -log log.water.pyscf.qmmm.plugin.1 -in in.water.pyscf.qmmm.plugin diff --git a/examples/QUANTUM/PySCF/data.water.pyscf b/examples/QUANTUM/PySCF/data.water.pyscf new file mode 100644 index 0000000000..01b2df64e6 --- /dev/null +++ b/examples/QUANTUM/PySCF/data.water.pyscf @@ -0,0 +1,51 @@ +LAMMPS data file via write_data, version 29 Sep 2021, timestep = 100000 + +6 atoms +2 atom types +4 bonds +1 bond types +2 angles +1 angle types + +0 31.35187685 xlo xhi +0 31.35187685 ylo yhi +0 31.35187685 zlo zhi + +Masses + +1 16 +2 1.008 + +Pair Coeffs # lj/cut/coul/long + +1 0.155425 3.16549 +2 0 0 + +Bond Coeffs # harmonic + +1 529.581 1.012 + +Angle Coeffs # harmonic + +1 37.95 113.24 + +Atoms + +1 1 1 -0.82 1.5026689324077283 5.863928543354579 3.792611836941212 +2 1 2 0.41 0.949608563235617 6.67792878850468 4.024506527262098 +3 1 2 0.41 2.542061320328983 6.034841727577386 3.822920239902413 +4 2 1 -0.82 7.984023446667837 11.359952374736203 6.435736949920884 +5 2 2 0.41 8.817227182546942 11.213723954966298 5.880386471510548 +6 2 2 0.41 7.933347339276622 12.376737335657605 6.56734173868257 + +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 diff --git a/examples/QUANTUM/PySCF/in.water.pyscf.qmmm b/examples/QUANTUM/PySCF/in.water.pyscf.qmmm new file mode 100644 index 0000000000..68a2c010e7 --- /dev/null +++ b/examples/QUANTUM/PySCF/in.water.pyscf.qmmm @@ -0,0 +1,46 @@ +# QMMM with PySCF - two water example + +units real +atom_style full + +pair_style lj/cut/coul/long 12 +pair_modify mix arithmetic + +bond_style harmonic +angle_style harmonic +dihedral_style none +improper_style none + +kspace_style pppm 1e-5 + +read_data data.water.pyscf + +# QM atoms are 1st water +# MM atoms are 2nd water + +group qm molecule 1 +group mm molecule 2 + +# remove bonds/angles between QM atoms +# set charges to zero on QM atoms + +delete_bonds qm multi remove special +set group qm charge 0.0 + +neighbor 2.0 bin +neigh_modify delay 0 every 1 check yes + +# QMMM dynamics + +timestep 2.0 + +fix 1 all nve + +fix 2 qm mdi/qmmm direct elements 8 1 +fix_modify 2 energy yes + +thermo_style custom step cpu temp ke evdwl ecoul epair emol elong & + f_2 pe etotal press + +thermo 1 +run 2 diff --git a/examples/QUANTUM/PySCF/in.water.pyscf.qmmm.plugin b/examples/QUANTUM/PySCF/in.water.pyscf.qmmm.plugin new file mode 100644 index 0000000000..4c9261d67b --- /dev/null +++ b/examples/QUANTUM/PySCF/in.water.pyscf.qmmm.plugin @@ -0,0 +1,49 @@ +# QMMM with PySCF +# adapted from pyscf_water/lmp.in for 2 water QMMM example + +units real +atom_style full + +pair_style lj/cut/coul/long 12 +pair_modify mix arithmetic + +bond_style harmonic +angle_style harmonic +dihedral_style none +improper_style none + +kspace_style pppm 1e-5 + +read_data data.water.pyscf + +# QM atoms are 1st water +# MM atoms are 2nd water + +group qm molecule 1 +group mm molecule 2 + +# remove bonds/angles between QM atoms +# set charges to zero on QM atoms + +delete_bonds qm multi remove special +set group qm charge 0.0 + +neighbor 2.0 bin +neigh_modify delay 0 every 1 check yes + +# QMMM dynamics + +timestep 2.0 + +fix 1 all nve + +fix 2 qm mdi/qmmm direct elements 8 1 +fix_modify 2 energy yes + +thermo_style custom step cpu temp ke evdwl ecoul epair emol elong & + f_2 pe etotal press + +thermo 1 + +mdi plugin pyscf_mdi mdi "-role ENGINE -name PySCF -method LINK" & + extra "-pbc yes" command "run 2" diff --git a/examples/QUANTUM/PySCF/pyscf_mdi.py b/examples/QUANTUM/PySCF/pyscf_mdi.py new file mode 100644 index 0000000000..459dc0d419 --- /dev/null +++ b/examples/QUANTUM/PySCF/pyscf_mdi.py @@ -0,0 +1,657 @@ +# MDI wrapper on PySCF quantum code + +# NOTE: Qs or issues to still address +# mm_radii is input to PySCF in Angstroms? +# add list of radii for all elements +# allow for changes in box size, e.g. every step for NPT +# PySCF can do DIRECT mode (LATTICE commands) for QMMM +# can it also do POTENTIAL mode (Coulomb potential at QM atoms) +# if so, add PySCF logic in evaluate() +# if not, remove support for MDI POTENTIAL_AT_NUCLEI command +# can PySCF return stress +# any other PySCF settings for users to set +# are the 3 below good default values +# is wiping out dm_previous sufficient to make it a new system +# should max_memory=10000 depend on # of MM or QM atoms +# doc what "mf" is in evaluate() +# add PySCF code for AIMD (no MM atoms) +# how to specify box for mixed BC, i.e. cell.dimension = 2 or 1 +# also need command-line options for those cases ? +# redirect PySCF output to a file ? + +import sys,time + +import numpy as np +from mpi4py import MPI +import MDI_Library as mdi + +from pyscf.gto import Mole +from pyscf.pbc.gto import Cell +from pyscf import qmmm +from pyscf.dft import RKS + +# -------------------------------------------- + +# atomic_number_to_radius converts atomic number to radius (Angstroms) +# Chemistry - A European Journal, (2009), 186-197, 15(1) + +atomic_number_to_radius = {1: 0.32, 6: 0.75, 7: 0.71, 8: 0.63, 17: 0.99} + +# ELEMENTS is from pyscf/pyscf/data/elements.py + +ELEMENTS = [ + 'H' , 'He', 'Li', 'Be', 'B' , 'C' , 'N' , 'O' , 'F' , 'Ne', + 'Na', 'Mg', 'Al', 'Si', 'P' , 'S' , 'Cl', 'Ar', 'K' , 'Ca', + 'Sc', 'Ti', 'V' , 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', + 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y' , 'Zr', + 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', + 'Sb', 'Te', 'I' , 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', + 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', + 'Lu', 'Hf', 'Ta', 'W' , 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', + 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', + 'Pa', 'U' , 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', + 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', + 'Rg', 'Cn', 'Nh', 'Fl', 'Mc', 'Lv', 'Ts', 'Og', +] + +# atomic_number_to_symbol converts atomic number to element symbol + +atomic_number_to_symbol = {} +for i,symbol in enumerate(ELEMENTS): + atomic_number_to_symbol[i+1] = symbol + +# -------------------------------------------- +# PySCF settings +# these are default values +# options() may override them +# -------------------------------------------- + +periodic = 1 +xcstr = "wb97x" +basis = "6-31+G**" + +# -------------------------------------------- +# persistent data +# -------------------------------------------- + +world = 0 +me = nprocs = 0 +exitflag = False + +AIMD = 0 +QMMM = 1 +mode = AIMD + +# QM inputs + +flag_qm_natoms = flag_mm_natoms = 0 +flag_box = flag_box_displ = 0 +flag_qm_elements = 0 +flag_qm_coords = flag_qm_potential = 0 +flag_mm_elements = 0 +flag_mm_coords = flag_mm_charges = 0 + +box = np.empty(9) +box_displ = np.empty(3) + +qm_natoms = 0 +qm_elements = None +qm_coords = None +qm_potential = None + +mm_natoms = 0 +mm_coords = None +mm_charges = None +mm_elements = None +mm_radii = None + +# QM outputs + +qm_pe = 0.0 +qm_stress = np.empty(9) +qm_forces = None +qm_charges = None + +mm_forces = None + +# PySCF internal data + +dm_previous_exists = 0 +dm_previous = None + +# -------------------------------------------- +# print error message and halt +# -------------------------------------------- + +def error(txt): + if me == 0: print("ERROR:",txt) + world.Abort() + +# -------------------------------------------- +# process non-MDI options to PySCF +# if this script is executed independently: +# args = command-line args +# if this script is invoked as a plugin library: +# args = passed via MDI +# -------------------------------------------- + +def options(args): + global periodic,xcstr,basis + + narg = len(args) + iarg = 0 + while iarg < narg: + if args[iarg] == "-pbc": + if iarg+1 > narg: error("Invalid PySCF command line args") + if args[iarg+1] == "yes": periodic = 1 + elif args[iarg+1] == "no": periodic = 0 + else: error("Invalid PySCF command line args") + iarg += 2 + elif args[iarg] == "-xcstr": + if iarg+1 > narg: error("Invalid PySCF command line args") + xcstr = args[iarg+1] + iarg += 2 + elif args[iarg] == "-basis": + if iarg+1 > narg: error("Invalid PySCF command line args") + basis = args[iarg+1] + iarg += 2 + else: error("Invalid PySCF command line args") + +# -------------------------------------------- +# operate as an engine +# -------------------------------------------- + +def mdi_engine(other_options): + global world + + # get the MPI intra-communicator for this engine + + world = mdi.MDI_MPI_get_world_comm() + me = world.Get_rank() + nprocs = world.Get_size() + + # process non-MDI command line args + + options(other_options) + + # confirm PySCF is being run as an engine + + role = mdi.MDI_Get_Role() + if not role == mdi.MDI_ENGINE: + error("Must run PySCF as an MDI engine") + + # supported MDI commands + + mdi.MDI_Register_Node("@DEFAULT") + mdi.MDI_Register_Command("@DEFAULT","EXIT") + + # driver --> engine + + mdi.MDI_Register_Command("@DEFAULT",">NATOMS") + mdi.MDI_Register_Command("@DEFAULT",">CELL") + mdi.MDI_Register_Command("@DEFAULT",">CELL_DISPL") + mdi.MDI_Register_Command("@DEFAULT",">ELEMENTS") + mdi.MDI_Register_Command("@DEFAULT",">COORDS") + mdi.MDI_Register_Command("@DEFAULT",">POTENTIAL_AT_NUCLEI") + mdi.MDI_Register_Command("@DEFAULT",">NLATTICE") + mdi.MDI_Register_Command("@DEFAULT",">LATTICE_ELEMENTS") + mdi.MDI_Register_Command("@DEFAULT",">CLATTICE") + mdi.MDI_Register_Command("@DEFAULT",">LATTICE") + + # engine --> driver + + mdi.MDI_Register_Command("@DEFAULT","LATTICE_FORCES") + mdi.MDI_Register_Command("@DEFAULT","NATOMS": + receive_qm_natoms(mdicomm) + + elif command == ">CELL": + receive_box(mdicomm) + + elif command == ">CELL_DISPL": + receive_box_displ(mdicomm) + + elif command == ">ELEMENTS": + receive_qm_elements(mdicomm) + + elif command == ">COORDS": + receive_qm_coords(mdicomm) + + elif command == ">POTENTIAL_AT_NUCLEI": + receive_qm_potential(mdicomm) + + elif command == ">NLATTICE": + receive_mm_natoms(mdicomm) + + elif command == ">LATTICE_ELEMENTS": + receive_mm_elements(mdicomm) + + elif command == ">CLATTICE": + receive_mm_coords(mdicomm) + + elif command == ">LATTICE": + receive_mm_charges(mdicomm) + + # MDI commands which retreive quantum results + # each may also trigger the quantum calculation + + elif command == "CELL data") + world.Bcast(box,root=0) + +# -------------------------------------------- +# receive simulation box displacement vector from driver +# -------------------------------------------- + +def receive_box_displ(mdicomm): + global flag_box_displ + flag_box_displ = 1 + + ierr = mdi.MDI_Recv(3,mdi.MDI_DOUBLE,mdicomm,buf=box_displ) + if ierr: error("MDI: >CELL_DISPL data") + world.Bcast(box_displ,root=0) + +# -------------------------------------------- +# receive QM atom coords from driver +# -------------------------------------------- + +def receive_qm_coords(mdicomm): + global flag_qm_coords + flag_qm_coords = 1 + + if not qm_natoms: error("Cannot MDI >COORDS if # of QM atoms = 0") + + ierr = mdi.MDI_Recv(3*qm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=qm_coords) + if ierr: error("MDI: >COORDS data") + world.Bcast(qm_coords,root=0) + +# -------------------------------------------- +# receive Coulomb potential at QM nuclei from driver +# -------------------------------------------- + +def receive_qm_potential(mdicomm): + global flag_qm_potential + flag_qm_potential = 1 + + if not qm_natoms: error("Cannot MDI >POTENTIAL_AT_NUCLEI if # of QM atoms = 0") + + ierr = mdi.MDI_Recv(qm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=qm_potential) + if ierr: error("MDI: >POTENTIAL_AT_NUCLEI data") + world.Bcast(qm_potential,root=0) + +# -------------------------------------------- +# receive QM atomic numbers from driver +# -------------------------------------------- + +def receive_qm_elements(mdicomm): + global flag_qm_elements + flag_qm_elements = 1 + + if not qm_natoms: error("Cannot MDI >ELEMENTS if # of QM atoms = 0") + + ierr = mdi.MDI_Recv(qm_natoms,mdi.MDI_INT,mdicomm,buf=qm_elements) + if ierr: error("MDI: >ELEMENTS data") + world.Bcast(qm_elements,root=0) + +# -------------------------------------------- +# receive count of MM atoms from driver +# -------------------------------------------- + +def receive_mm_natoms(mdicomm): + global flag_mm_natoms,mm_natoms + flag_mm_natoms = 1 + + mm_natoms = mdi.MDI_Recv(1,mdi.MDI_INT,mdicomm) + mm_natoms = world.bcast(mm_natoms,root=0) + allocate("mm") + +# -------------------------------------------- +# receive MM atomic numbers from driver +# -------------------------------------------- + +def receive_mm_elements(mdicomm): + global flag_mm_elements + flag_mm_elements = 1 + + if not mm_natoms: error("Cannot MDI >LATTICE_ELEMENTS if # of MM atoms = 0") + + ierr = mdi.MDI_Recv(mm_natoms,mdi.MDI_INT,mdicomm,buf=mm_elements) + if ierr: error("MDI: >LATTICE_ELEMENTS data") + world.Bcast(mm_elements,root=0) + +# -------------------------------------------- +# receive MM atom coords from driver +# -------------------------------------------- + +def receive_mm_coords(mdicomm): + global flag_mm_coords + flag_mm_coords = 1 + + if not mm_natoms: error("Cannot MDI >CLATTICE if # of MM atoms = 0") + + ierr = mdi.MDI_Recv(3*mm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=mm_coords) + if ierr: error("MDI: >CLATTICE data") + world.Bcast(mm_coords,root=0) + +# -------------------------------------------- +# receive charge on MM atoms from driver +# -------------------------------------------- + +def receive_mm_charges(mdicomm): + global flag_mm_charges + flag_mm_charges = 1 + + if not mm_natoms: error("Cannot MDI >LATTICE if # of MM atoms = 0") + + ierr = mdi.MDI_Recv(mm_natoms,mdi.MDI_DOUBLE,mdicomm,buf=mm_charges) + if ierr: error("MDI: >LATTICE data") + world.Bcast(mm_charges,root=0) + +# -------------------------------------------- +# allocate persistent data for QM or MM atoms +# called when qm_natoms or mm_natoms is reset by MDI driver +# -------------------------------------------- + +def allocate(which): + global qm_elements,qm_coords,qm_potential,qm_forces,qm_charges + global mm_elements,mm_coords,mm_charges,mm_forces + + if which == "qm": + n = qm_natoms + qm_elements = np.empty(n,dtype=np.int32) + qm_coords = np.empty((n,3)) + qm_potential = np.empty(n) + qm_forces = np.empty((n,3)) + qm_charges = np.empty(n) + + if which == "mm": + n = mm_natoms + mm_elements = np.empty(n,dtype=np.int32) + mm_coords = np.empty((n,3)) + mm_charges = np.empty(n) + mm_forces = np.empty((n,3)) + +# -------------------------------------------- +# perform a quantum calculation via PySCF +# -------------------------------------------- + +def evaluate(): + global mode; + global flag_qm_natoms,flag_mm_natoms + global flag_box,flag_box_displ + global flag_qm_elements,flag_qm_coords,flag_qm_potential + global flag_mm_elements,flag_mm_coords,flag_mm_charges + global qm_pe,qm_stress,qm_forces,qm_charges + global mm_forces + global dm_previous_exists,dm_previous + + # just return if the QM system was already evaluated + # happens when multiple results are requested by driver + + any_flag = flag_qm_natoms + flag_mm_natoms + flag_box + flag_box_displ + \ + flag_qm_elements + flag_qm_coords + flag_qm_potential + \ + flag_mm_elements + flag_mm_coords + flag_mm_charges + if not any_flag: return + + # if any of these MDI commands received from LAMMPS, + # treat it as a brand new system + + new_system = 0 + if flag_qm_natoms or flag_mm_natoms: new_system = 1 + if flag_qm_elements or flag_mm_elements: new_system = 1 + if new_system: + if flag_mm_natoms or flag_qm_potential: mode = QMMM + else: mode = AIMD + dm_previous_exists = 0 + + # if new system, error check that all needed MDI calls have been made + + if new_system: + if periodic and not flag_box: + error("Simulation box not specified for periodic system") + if not flag_qm_natoms: error("QM atom count not specified") + if not flag_qm_elements or not flag_qm_coords: + error("QM atom properties not fully specified") + if flag_mm_natoms: + if not flag_mm_elements or not flag_mm_coords or not flag_mm_charges: + error("MM atom properties not fully specified") + + # PySCF inputs for QM and MM atoms + # box, qm_coords, mm_coords must be converted to Angstroms + + bohr_to_angstrom = mdi.MDI_Conversion_factor("bohr","angstrom") + + box_A = box * bohr_to_angstrom + qm_coords_A = qm_coords * bohr_to_angstrom + mm_coords_A = mm_coords * bohr_to_angstrom + + qm_symbols = [atomic_number_to_symbol[anum] for anum in qm_elements] + mm_radii = [atomic_number_to_radius[anum] for anum in mm_elements] + + #pos2str = lambda pos: " ".join([str(x) for x in qm_coords_A]) + #atom_str = [f"{a} {pos2str(pos)}\n" for a,pos in zip(qm_symbols,qm_coords_A)] + + clines = ["%s %20.16g %20.16g %20.16g" % (symbol,xyz[0],xyz[1],xyz[2]) + for symbol,xyz in zip(qm_symbols,qm_coords_A)] + atom_str = "\n".join(clines) + + if periodic: + edge_vec = "%20.16g %20.16g %20.16g" + box_str = "%s\n%s\n%s" % (edge_vec,edge_vec,edge_vec) + box_str = box_str % \ + (box_A[0],box_A[1],box_A[2],box_A[3],box_A[4],box_A[5],box_A[6],box_A[7],box_A[8]) + + print("ATOM STR:",atom_str) + print("BOX STR:",box_str) + print("MM COORDS:",mm_coords) + print("MM CHARGES:",mm_charges) + print("MM RADII:",mm_radii) + + # build PySCF system + # use Cell for periodic, Mole for non-periodic + + if periodic: + cell = Cell() + cell.atom = atom_str + cell.a = box_str + cell.max_memory = 10000 + cell.basis = basis + cell.build() + else: + mol = Mole() + mol.atom = atom_str + mol.max_memory = 10000 + mol.basis = basis + mol.build() + + # QMMM with QM and MM atoms + # mf = ??? + # qm_pe = QM energy + QM/MM energy + # QM energy = QM_nuclear/QM_nuclear + electron/QM_nuclear + electron/electron + # QM/MM energy = QM_nuclear/MM_charges + electron/MM_charges + # qm_forces = QM forces = same 3 terms + # mm_forces = QM/MM forces = same 2 terms + # dm = molecular orbitals (wave functions) for system + + if mode == QMMM: + if periodic: mf = RKS(cell,xc=xcstr) + else: mf = RKS(mol,xc=xcstr) + mf = qmmm.mm_charge(mf,mm_coords,mm_charges,mm_radii) + + if dm_previous_exists: + qm_pe = mf.kernel(dm0=dm_previous) + else: + qm_pe = mf.kernel() + + mf_grad = mf.nuc_grad_method() + qm_forces = -mf_grad.kernel() + dm = mf.make_rdm1() + mm_forces = -(mf_grad.grad_nuc_mm() + mf_grad.contract_hcore_mm(dm)) + + dm_previous_exists = 1 + dm_previous = dm + + # AIMD with only QM atoms + + elif mode == AIMD: + pass + + # clear flags for all MDI commands for next QM evaluation + + flag_qm_natoms = flag_mm_natoms = 0 + flag_box = flag_box_displ = 0 + flag_qm_elements = 0 + flag_qm_coords = flag_qm_potential = 0 + flag_mm_elements = 0 + flag_mm_coords = flag_mm_charges = 0 + +# -------------------------------------------- +# function called by MDI driver +# only when it invokes pyscf_mdi.py as a plugin +# -------------------------------------------- + +def MDI_Plugin_init_pyscf_mdi(plugin_state): + + # other_options = all non-MDI args + # -mdi arg is processed and stripped internally by MDI + + print("PLUGIN entry") + + other_options = [] + + mdi.MDI_Set_plugin_state(plugin_state) + narg = mdi.MDI_Plugin_get_argc() + + for iarg in range(narg): + arg = mdi.MDI_Plugin_get_arg(iarg) + other_options.append(arg) + + # start running as an MDI engine + + mdi_engine(other_options) + +# -------------------------------------------- +# main program +# invoked when MDI driver and pyscf_mdi.py +# are run as independent executables +# -------------------------------------------- + +if __name__== "__main__": + + # mdi_option = single arg in quotes that follows -mdi + # other_options = all non-MDI args + + mdi_option = "" + other_options = [] + + narg = len(sys.argv) + args = sys.argv + + iarg = 1 + while iarg < narg: + arg = args[iarg] + if arg == "-mdi" or arg == "--mdi": + if narg > iarg+1: mdi_option = sys.argv[iarg+1] + else: error("PySCF -mdi argument not provided") + iarg += 1 + else: other_options.append(arg) + iarg += 1 + + if not mdi_option: error("PySCF -mdi option not provided") + + # call MDI_Init with just -mdi option + + mdi.MDI_Init(mdi_option) + + # start running as an MDI engine + + mdi_engine(other_options) diff --git a/examples/QM/README b/examples/QUANTUM/README similarity index 76% rename from examples/QM/README rename to examples/QUANTUM/README index 6028d2be8d..009ae3ffc1 100644 --- a/examples/QM/README +++ b/examples/QUANTUM/README @@ -5,7 +5,16 @@ LATTE = semi-empirical tight-binding code from LANL https://www.osti.gov/biblio/1526907-los-alamos-transferable-tight-binding-energetics-latte-version https://github.com/lanl/LATTE -To be added later (as of Aug 2022): +PySCF = ??? from Caltech + add link + +----------------------------------------------------- + +To be added later (as of Jan 2023): + +NWChem = computational chemistry code from PNNL + focus here is on DFT portion of NWChem = PWDFT + https://www.nwchem-sw.org Quantum Espresso (QE) = DFT code for materials modeling https://www.quantum-espresso.org/ @@ -16,6 +25,3 @@ DFT-FE = real-space DFT code from U Michigan INQ = DFT code from LLNL https://github.com/LLNL/inq -NWChem = computational chemistry code from PNNL - focus here is on DFT portion of NWChem - https://www.nwchem-sw.org