diff --git a/examples/QUANTUM/NWChem/README b/examples/QUANTUM/NWChem/README new file mode 100644 index 0000000000..3d570cb609 --- /dev/null +++ b/examples/QUANTUM/NWChem/README @@ -0,0 +1,119 @@ +# Test runs of QMMM with LAMMPS and NWChem + +Step 1: build LAMMPS +Step 2: download/build the MDI code coupling package +Step 3: download/build NWChem PWDFT +Step 4: run 2-water QMMM problem for a few steps + +--------------------------------- +--------------------------------- + +Step 1: build LAMMPS + +The molecule package is needed for the 2-water test +problem. Copy the final LAMMPS executable into the +examples/QUANTUM/NWChem directory. + +Traditional make: + +% cd ~/lammps/lib/mdi +% python Install.py -m mpi +% cd ~/lammps/src +% make yes-mdi yes-molecule +% make -j mpi +% cp lmp_mpi ~/lammps/examples/QUANTUM/NWChem + +CMake: + +% cd ~/lammps +% mkdir build; cd build +% cmake -D PKG_MDI=yes -D PKG_MOLECULE=yes ../cmake +% make -j +% cp lmp ~/lammps/examples/QUANTUM/NWChem/lmp_mpi + +--------------------------------- +--------------------------------- + +Step 2: install the MDI code coupling package + +(a) grab the MDI Git repo + +% mkdir mdi; cd mdi +% git clone git@github.com:MolSSI-MDI/MDI_Library.git git + +(b) build MDI + +% cd mdi/git +% mkdir build; cd build +% cmake .. +% make -j + +(c) Add something similar to the following to your .bashrc or .cshrc +file so that Python can find MDI: + +For bash: + +% export PYTHONPATH="$PYTHONPATH:/home/sjplimp/mdi/git" +% hash -r + +For (t)csh: + +% setenv PYTHONPATH ${PYTHONPATH}:/home/sjplimp/mdi/git +% rehash + +(d) Check that you can import the 3 Python modules which the script +that wraps PySCF will need: + +% python +>>> import numpy as np +>>> from mpi4py import MPI +>>> import MDI_Library as mdi + +--------------------------------- +--------------------------------- + +Step 3: download/build NWChem PWDFT + +(a) grab the PWDFT Git repo + +% mkdir nwchem; cd nwchem +% git clone git@github.com:ebylaska/PWDFT.git PWDFT + +(b) build PWDFT + +% cd nwchem/PWDFT +% cd build_library; rm -r * +% cmake ../Nwpw -DMAKE_LIBRARY=true -DCMAKE_POSITION_INDEPENDENT_CODE=ON +% make -j # should produce libpwdft.so in build_library + +(c) Add something similar to the following to your .bashrc or .cshrc +file so that the Python wrapper script can find libpwdft.so: + +For bash: + +% export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/home/sjplimp/nwchem/PWDFT/build_library" +% hash -r + +For (t)csh: + +% setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:/home/sjplimp/nwchem/PWDFT/build_library +% rehash + +--------------------------------- +--------------------------------- + +Step 4: run the 2-water QMMM problem for a few steps + +% cd ~/lammps/examples/QUANTUM/NWChem + +lmp_mpi -mdi "-name LMP -role DRIVER -method TCP -port 8021" -log log.water.nwchem.qmmm.tcp.1 -in in.water.nwchem.qmmm & + +python nwchem_mdi.py -mdi "-name NWChem -role ENGINE -method TCP -port 8021 -hostname localhost" template.water.nw water.dimer.nw log.water.pwdft.qmmm.tcp.1 + +# Run with MPI: 1 proc each + +mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" -log log.water.nwchem.qmmm.mpi.1 -in in.water.nwchem.qmmm : -np 1 python nwchem_mdi.py -mdi "-name NWChem -role ENGINE -method MPI" template.water.nw water.dimer.nw log.water.pwdft.qmmm.mpi.1 + +# Run in plugin mode: 1 proc + +lmp_mpi -mdi "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/work_qm" -log log.water.nwchem.qmmm.plugin.1 -in in.water.nwchem.qmmm.plugin diff --git a/examples/QUANTUM/NWChem/data.water.nwchem b/examples/QUANTUM/NWChem/data.water.nwchem new file mode 100644 index 0000000000..cd21c69f52 --- /dev/null +++ b/examples/QUANTUM/NWChem/data.water.nwchem @@ -0,0 +1,48 @@ +LAMMPS data file for water dimer in large box - for QMMM with different types + +6 atoms +4 atom types +4 bonds +1 bond types +2 angles +1 angle types + +-6.879301 6.879301 xlo xhi +-6.879301 6.879301 ylo yhi +-6.879301 6.879301 zlo zhi + +Masses + +1 15.99491 +2 1.008 +3 15.99491 +4 1.008 + +Bond Coeffs + +1 554.25 1.0 + +Angle Coeffs + +1 47.744 109.4 + +Atoms + +1 1 1 -0.8476 0.161560 -0.052912 0.033173 +2 1 2 0.4238 0.803054 0.369132 -0.511660 +3 1 2 0.4238 -0.325571 -0.669574 -0.488560 +4 2 3 -0.8476 0.021259 0.506771 2.831278 +5 2 4 0.4238 -0.721039 1.083100 2.758378 +6 2 4 0.4238 0.158220 0.181883 1.945696 + +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/NWChem/in.water.nwchem.qmmm b/examples/QUANTUM/NWChem/in.water.nwchem.qmmm new file mode 100644 index 0000000000..4b44002f58 --- /dev/null +++ b/examples/QUANTUM/NWChem/in.water.nwchem.qmmm @@ -0,0 +1,70 @@ +# QMMM with NWChem + +units real +atom_style full + +bond_style harmonic +angle_style harmonic + +read_data data.water.nwchem + +# QM atoms are 1st water +# MM atoms are 2nd water + +group qm molecule 1 +group mm molecule 2 + +# remove bonds/angles between QM atoms + +delete_bonds qm multi remove special + +# pair style must define stand-alone short-range Coulombics +# must specify mixing explicitly b/c hybrid/overlay +# QM O,H = types 1,2 +# MM O,H = types 3,4 +# QM O,H atoms do not LJ interact with each other +# only MM O atoms LJ interact with other b/c MM H is zero +# MM/QM O do LJ interact with each other, same as pair of MM O atoms +# MM O and QM H do LJ interact with each other with non-zero H epsilon = 0.044 +# geometric mixing for epsilon, arithmetic for sigma +# this is to provide stability for QM H atoms + +# mixing only for MM-O/QM-O and MM-O/QM-H + +pair_style hybrid/overlay lj/cut 6.0 coul/cut 6.0 +pair_coeff 1 1 lj/cut 0.0 3.165558 +pair_coeff 2 2 lj/cut 0.0 0.7 +pair_coeff 3 3 lj/cut 0.155394 3.165558 +pair_coeff 4 4 lj/cut 0.0 0.7 +pair_coeff 1 3 lj/cut 0.155394 3.165558 +pair_coeff 2 3 lj/cut 0.08268818537130924 1.932779 +pair_coeff * * coul/cut + +neighbor 1.0 bin +neigh_modify delay 0 every 1 check yes + +# QMMM dynamics + +timestep 0.1 + +fix 1 all nve + +fix 2 qm mdi/qmmm potential elements 8 1 8 1 +#fix 2 qm nwchem template.water.nw water.dimer2.nw & +# log.pwdft.water.dimer O H O H +fix_modify 2 energy yes + +thermo_style custom step cpu temp ke evdwl ecoul epair emol elong & + f_2 pe etotal press + +# convert dump file forces to Hartree/Bohr for comparison to NWChem + +variable fx atom fx/1185.8 +variable fy atom fy/1185.8 +variable fz atom fz/1185.8 + +dump 1 all custom 1 dump.water.nwchem.qmmm id x y z q v_fx v_fy v_fz +dump_modify 1 sort id format float "%20.16g" + +thermo 1 +run 2 diff --git a/examples/QUANTUM/NWChem/in.water.nwchem.qmmm.plugin b/examples/QUANTUM/NWChem/in.water.nwchem.qmmm.plugin new file mode 100644 index 0000000000..4460a04dc8 --- /dev/null +++ b/examples/QUANTUM/NWChem/in.water.nwchem.qmmm.plugin @@ -0,0 +1,72 @@ +# QMMM with NWChem + +units real +atom_style full + +bond_style harmonic +angle_style harmonic + +read_data data.water.nwchem + +# QM atoms are 1st water +# MM atoms are 2nd water + +group qm molecule 1 +group mm molecule 2 + +# remove bonds/angles between QM atoms + +delete_bonds qm multi remove special + +# pair style must define stand-alone short-range Coulombics +# must specify mixing explicitly b/c hybrid/overlay +# QM O,H = types 1,2 +# MM O,H = types 3,4 +# QM O,H atoms do not LJ interact with each other +# only MM O atoms LJ interact with other b/c MM H is zero +# MM/QM O do LJ interact with each other, same as pair of MM O atoms +# MM O and QM H do LJ interact with each other with non-zero H epsilon = 0.044 +# geometric mixing for epsilon, arithmetic for sigma +# this is to provide stability for QM H atoms + +# mixing only for MM-O/QM-O and MM-O/QM-H + +pair_style hybrid/overlay lj/cut 6.0 coul/cut 6.0 +pair_coeff 1 1 lj/cut 0.0 3.165558 +pair_coeff 2 2 lj/cut 0.0 0.7 +pair_coeff 3 3 lj/cut 0.155394 3.165558 +pair_coeff 4 4 lj/cut 0.0 0.7 +pair_coeff 1 3 lj/cut 0.155394 3.165558 +pair_coeff 2 3 lj/cut 0.08268818537130924 1.932779 +pair_coeff * * coul/cut + +neighbor 1.0 bin +neigh_modify delay 0 every 1 check yes + +# QMMM dynamics + +timestep 0.1 + +fix 1 all nve + +fix 2 qm mdi/qmmm potential elements 8 1 8 1 +fix_modify 2 energy yes + +thermo_style custom step cpu temp ke evdwl ecoul epair emol elong & + f_2 pe etotal press + +# convert dump file forces to Hartree/Bohr for comparison to NWChem + +variable fx atom fx/1185.8 +variable fy atom fy/1185.8 +variable fz atom fz/1185.8 + +dump 1 all custom 1 dump.water.dimer.qmmm.plugin & + id x y z q v_fx v_fy v_fz +dump_modify 1 sort id format float "%20.16g" + +thermo 1 + +mdi plugin nwchem_mdi mdi "-role ENGINE -name NWChem -method LINK" & + extra "template.water.nw water.dimer.nw log.water.pwdft.qmmm.plugin.1" & + command "run 2" diff --git a/examples/QUANTUM/NWChem/nwchem_mdi.py b/examples/QUANTUM/NWChem/nwchem_mdi.py new file mode 100644 index 0000000000..0208a269a4 --- /dev/null +++ b/examples/QUANTUM/NWChem/nwchem_mdi.py @@ -0,0 +1,662 @@ +# MDI wrapper on NWChem PWDFT code + +# NOTE: Qs or issues to still address +# test if works for both AIMD and QMMM +# can series of problem be run via lib interface input_filename() ? +# how does PBC vs non-PBC work, just box size in NWC input file +# or maybe other settings in that file? +# can NWChem return stress? +# can NWChem do DIRECT mode? +# any options for 2d or 1d periodic? +# allow for box size changes, e.g. every step for NPT +# check NWC func call error returns ? + +import sys,time +from ctypes import * + +import numpy as np +from mpi4py import MPI +import MDI_Library as mdi + +# -------------------------------------------- + +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 + +# -------------------------------------------- +# global data +# -------------------------------------------- + +world = 0 +me = nprocs = 0 + +if MPI._sizeof(MPI.Comm) == sizeof(c_int): + MPI_Comm = c_int +else: + MPI_Comm = c_void_p + +exitflag = False + +AIMD = 0 +QMMM = 1 +mode = AIMD + +# NWChem PWDFT library + +libname = "libpwdft.so" +libpwdft = None + +# QM inputs + +nw_template = "" +nw_infile = "" +nw_outfile = "" + +flag_qm_natoms = flag_mm_natoms = 0 +flag_box = flag_box_displs = 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 + +# QM outputs + +qm_pe = 0.0 +qm_stress = np.empty(9) +qm_forces = None +qm_charges = None + +mm_forces = None + +# -------------------------------------------- +# print error message and halt +# -------------------------------------------- + +def error(txt): + if me == 0: print("ERROR:",txt) + world.Abort() + +# -------------------------------------------- +# process non-MDI options to PWDFT +# 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(other_options): + global nw_template,nw_infile,nw_outfile + if len(other_options) != 3: + error("Invalid args to NWChem wrapper: template_file infile outfile") + nw_template = other_options[0] + nw_infile = other_options[1] + nw_outfile = other_options[2] + +# -------------------------------------------- +# operate as an engine +# -------------------------------------------- + +def mdi_engine(other_options): + global world,MPI_Comm,libpwdft + + # 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 PWDFT is being run as an engine + + role = mdi.MDI_Get_Role() + if not role == mdi.MDI_ENGINE: + error("Must run NWChem 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_DISPLS") + 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 NWChem PWDFT +# -------------------------------------------- + +def evaluate(): + global mode + global flag_qm_natoms,flag_mm_natoms + global flag_box,flag_box_displs + 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 + + # 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_displs + \ + 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 + + # if new system, error check that all needed MDI calls have been made + + if new_system: + if not flag_qm_natoms: error("QM atom count not specified") + if not flag_qm_elements or not flag_qm_coords or not flag_qm_potential: + error("QM atom properties not fully specified") + + # setup new system within PWDFT + + if new_system: pwdft_initialize() + + # QMMM with QM and MM atoms + + world_ptr = MPI._addressof(world) + c_world = MPI_Comm.from_address(world_ptr) + c_qm_pe = c_double(qm_pe) + c_nw_outfile = nw_outfile.encode() + + if mode == QMMM: + print("QMMM minimizer") + time1 = time.time() + nwerr = libpwdft.\ + c_lammps_pspw_qmmm_minimizer_filename(c_world,qm_coords,qm_potential, + qm_forces,qm_charges,byref(c_qm_pe), + False,True,c_nw_outfile) + # NOTE: check nwerr return? + qm_pe = c_qm_pe.value + time2 = time.time() + print("DONE QMMM minimizer",nwerr,time2-time1) + print("PE",qm_pe) + print("FORCE",qm_forces) + print("CHARGES",qm_charges) + + # AIMD with only QM atoms + + elif mode == AIMD: + print("AIMD minimizer") + time1 = time.time() + nwerr = libpwdft.\ + c_lammps_pspw_aimd_minimizer_filename(c_world,qm_coords,qm_forces, + byref(c_qm_pe),c_nw_outfile) + # NOTE: check nwerr return? + qm_pe = c_qm_pe.value + time2 = time.time() + print("DONE AIMD minimizer",nwerr,time2-time1) + + # clear flags for all MDI commands for next QM evaluation + + flag_qm_natoms = flag_mm_natoms = 0 + flag_box = flag_box_displs = 0 + flag_qm_elements = 0 + flag_qm_coords = flag_qm_potential = 0 + flag_mm_elements = 0 + flag_mm_coords = flag_mm_charges = 0 + +# -------------------------------------------- +# load PWDFT lib and set ctypes signatures for function calls +# set ctypes signatures for 3 function calls to PWDFT lib +# -------------------------------------------- + +def pwdft_load(): + global libpwdft + + libpwdft = CDLL(libname,RTLD_GLOBAL) + + libpwdft.c_lammps_pspw_input_filename.restype = None + libpwdft.c_lammps_pspw_input_filename.argtypes = \ + [MPI_Comm, c_char_p, c_char_p] + + nparray = np.ctypeslib.ndpointer(dtype=np.float64,ndim=2,flags="C_CONTIGUOUS") + npvector = np.ctypeslib.ndpointer(dtype=np.float64,ndim=1,flags="C_CONTIGUOUS") + + libpwdft.c_lammps_pspw_qmmm_minimizer_filename.restype = c_int + libpwdft.c_lammps_pspw_qmmm_minimizer_filename.argtypes = \ + [MPI_Comm, nparray, npvector, nparray, npvector, POINTER(c_double), + c_bool, c_bool, c_char_p] + + libpwdft.c_lammps_pspw_aimd_minimizer_filename.restype = c_int + libpwdft.c_lammps_pspw_aimd_minimizer_filename.argtypes = \ + [MPI_Comm, nparray, nparray, POINTER(c_double), c_char_p] + +# -------------------------------------------- +# create PWDFT input file with box and list of atoms +# invoke PWDFT function to read it +# -------------------------------------------- + +def pwdft_initialize(): + + # box, qm_coords, mm_coords must be converted to Angstroms + + angstrom_to_bohr = mdi.MDI_Conversion_factor("angstrom","bohr") + bohr_to_angstrom = 1.0 / angstrom_to_bohr + + box_A = box * bohr_to_angstrom + qm_coords_A = qm_coords * bohr_to_angstrom + + # proc 0 reads template file, writes PWDFT input file + + if me == 0: + lines = open(nw_template,'r').readlines() + + fp = open(nw_infile,'w') + + for line in lines: + word = line.strip() + if word == "GEOMINSERT": + print("geometry noautosym noautoz nocenter",file=fp); + print("system crystal cartesian",file=fp) + print("lattice_vectors",file=fp) + print("%20.16g %20.16g %20.16g" % (box_A[0],box_A[1],box_A[2]),file=fp) + print("%20.16g %20.16g %20.16g" % (box_A[3],box_A[4],box_A[5]),file=fp) + print("%20.16g %20.16g %20.16g" % (box_A[6],box_A[7],box_A[8]),file=fp) + print("end\n",file=fp) + + for i in range(qm_natoms): + symbol = atomic_number_to_symbol[qm_elements[i]] + print("%s %20.16g %20.16g %20.16g" % + (symbol,qm_coords_A[i][0],qm_coords_A[i][1],qm_coords_A[i][2]), + file=fp) + print("end\n",file=fp) + + else: print(line,file=fp,end="") + + fp.close() + + # all procs call pspw_input_filename() which processes input file + # performs initial QM calculation within PWDFT + + world_ptr = MPI._addressof(world) + c_world = MPI_Comm.from_address(world_ptr) + infile = nw_infile.encode() + outfile = nw_outfile.encode() + + print("INPUT filename") + time1 = time.time() + nwerr = libpwdft.c_lammps_pspw_input_filename(c_world,infile,outfile) + time2 = time.time() + print("DONE INPUT filename",nwerr,time2-time1) + +# -------------------------------------------- +# function called by MDI driver +# only when it invokes pyscf_mdi.py as a plugin +# -------------------------------------------- + +def MDI_Plugin_init_nwchem_mdi(plugin_state): + + # other_options = all non-MDI args + # -mdi arg is processed and stripped internally by MDI + + 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("NWChem -mdi argument not provided") + iarg += 1 + else: other_options.append(arg) + iarg += 1 + + if not mdi_option: error("NWChem -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/QUANTUM/NWChem/template.water.nw b/examples/QUANTUM/NWChem/template.water.nw new file mode 100644 index 0000000000..5cd863593b --- /dev/null +++ b/examples/QUANTUM/NWChem/template.water.nw @@ -0,0 +1,17 @@ +Title "LAMMPS wrapping of PWDFT" + +memory 1900 mb + +echo + +GEOMINSERT + +nwpw + xc pbe + cutoff 30.0 + 2d-hcurve + tolerances 1.0e-9 1.0-9 + apc on +end + +task pspw gradient diff --git a/examples/QUANTUM/README b/examples/QUANTUM/README index 009ae3ffc1..19925b30ac 100644 --- a/examples/QUANTUM/README +++ b/examples/QUANTUM/README @@ -1,5 +1,5 @@ -Each of the directories shows how to use LAMMPS in tandem with a -specific quantum code +Each directory shows how to use LAMMPS in tandem with a specific +quantum code: LATTE = semi-empirical tight-binding code from LANL https://www.osti.gov/biblio/1526907-los-alamos-transferable-tight-binding-energetics-latte-version @@ -8,14 +8,14 @@ LATTE = semi-empirical tight-binding code from LANL PySCF = ??? from Caltech add link +NWChem = computational chemistry code from PNNL + focus here is on DFT portion of NWChem = PWDFT library + https://www.nwchem-sw.org + ----------------------------------------------------- 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/