Qmmm

examples/QMMM/2E4E_RHF-DFT-QMMM_energy.inp

@@ -0,0 +1,19 @@
+[input]
+system=2E4E.pdb 1 2 5 6 7 8 9 3 4 10 22 11 23
+charge=1
+runtype=energy
+functional=bhhlyp
+basis=6-31gs
+method=hf
+qmmm_flag=True
+
+[guess]
+type=huckel
+
+[scf]
+multiplicity=1
+type=rhf
+
+[dftgrid]
+rad_type=becke
+

examples/QMMM/ala.inp

@@ -0,0 +1,19 @@
+[input]
+system=ala.pdb 9 10 17 18 19
+charge=0
+runtype=energy
+functional=bhhlyp
+basis=6-31gs
+method=hf
+qmmm_flag=True
+
+[guess]
+type=huckel
+
+[scf]
+multiplicity=1
+type=rhf
+
+[dftgrid]
+rad_type=becke
+

examples/QMMM/ala.pdb

@@ -0,0 +1,23 @@
+HEADER    ala
+COMPND
+SOURCE
+ATOM      1  CH3 ACE     1      -0.028   0.101  -0.091
+ATOM      2  C   ACE     1      -0.023   0.157   1.419
+ATOM      3  O   ACE     1       1.014   0.404   2.021
+ATOM      4  H1  ACE     1       0.978   0.305  -0.455
+ATOM      5  H2  ACE     1      -0.715   0.851  -0.477
+ATOM      6  H3  ACE     1      -0.339  -0.891  -0.414
+ATOM      7  N   ALA     2      -1.183  -0.077   2.026
+ATOM      8  CA  ALA     2      -1.370  -0.108   3.477
+ATOM      9  C   ALA     2      -2.450  -1.132   3.872
+ATOM     10  O   ALA     2      -3.296  -1.508   3.068
+ATOM     11  H   ALA     2      -1.986  -0.315   1.463
+ATOM     12  HA  ALA     2      -0.433  -0.408   3.951
+ATOM     13  CB  ALA     2      -1.730   1.305   3.959
+ATOM     14  HB1 ALA     2      -1.842   1.311   5.044
+ATOM     15  HB2 ALA     2      -2.665   1.629   3.500
+ATOM     16  HB3 ALA     2      -0.935   2.000   3.685
+ATOM     17  N   NH2     3      -2.452  -1.588   5.119
+ATOM     18  H1  NH2     3      -1.760  -1.272   5.777
+ATOM     19  H2  NH2     3      -3.162  -2.253   5.377
+END

include/oqp.h

@@ -77,7 +77,7 @@ struct dft_parameters {
     bool dft_wt_der;
 };
 
- struct tddft_parameters {
+struct tddft_parameters {
     int64_t nstate;
     int64_t target_state;
     int64_t maxvec;
@@ -123,6 +123,7 @@ struct control_parameters {
     double    esp_constr;
     bool      basis_set_issue;
     double    conf_print_threshold;
+    bool     qmmm_flag;
 };
 
 struct mpi_communicator {

pyoqp/oqp/library/ints_1e.py

@@ -1,5 +1,11 @@
 import oqp
+import numpy as np
+from oqp.utils.qmmm import openmm_potential,openmm_energy
 
 def ints_1e(mol):
     """Compute a set of one-electron integrals"""
+    if mol.config['input']['qmmm_flag']:
+       current_xyz = mol.get_system().reshape((-1, 3))
+       mol.data["OQP::mm_potential"]=np.array(openmm_potential(current_xyz)).tolist()
+       mol.data["OQP::mm_energy"]=openmm_energy()
     oqp.int1e(mol)

pyoqp/oqp/library/libdlfind.py

@@ -8,6 +8,7 @@
 import numpy as np
 from oqp.library.libscipy import Optimizer
 from oqp.utils.file_utils import dump_data, dump_log
+import oqp.utils.qmmm as qmmm
 from libdlfind import dl_find
 from libdlfind.callback import (
     dlf_get_gradient_wrapper,
@@ -270,3 +271,98 @@ def check_convergence(self):
                 dump_log(self.mol,
                          title='PyOQP: Geometry Optimization Has Not Converged. Reached The Maximum Iteration')
                 raise StopIteration
+
+class DLFindQMMM(DLFinder):
+    """
+    OQP DL-Find state specific optimization class
+
+    """
+    def __init__(self, mol):
+        super().__init__(mol)
+
+        dump_log(self.mol, title='PyOQP: DL-FIND Local Minimum QM/MM Optimization', section='dlf')
+
+        self.ref_coord,self.mass,self.natom,self.atoms = qmmm.openmm_info()
+        self.pre_coord = np.array(self.ref_coord).reshape((-1))
+
+        self.dlf_get_params = make_dlf_get_params(
+            coords=self.pre_coord.tolist(),
+            printl=self.printl,
+            icoord=self.icoord,
+            iopt=self.iopt,
+            imultistate=self.ims,
+            maxcycle=self.maxit,
+            tolerance=1e-20,
+            tolerance_e=1e-20,
+        )
+
+
+    def one_step(self, coordinates):
+        # flatten coord
+        coordinates = coordinates.reshape(-1)
+
+        # add iteration
+        self.itr += 1
+
+        dump_log(self.mol, title='PyOQP: QM/MM Geometry Optimization Step %s' % self.itr)
+
+        # update coordinates
+        qmmm.openmm_update_system(coordinates)
+        current_xyz=coordinates.reshape((-1,3))
+#       coordinates_qm = np.array (())
+#       for i in qmmm.qm_atoms:
+#           coordinates_qm=np.append(coordinates_qm,current_xyz[i][0])
+#           coordinates_qm=np.append(coordinates_qm,current_xyz[i][1])
+#           coordinates_qm=np.append(coordinates_qm,current_xyz[i][2])
+
+        num_atoms, x, y, z, q, mass = qmmm.openmm_system()
+        coordinates_qm=np.array(x + y + z).reshape((3, num_atoms)).T.reshape(-1)
+
+        self.mol.update_system(coordinates_qm)
+
+        if self.itr == 1:
+            do_init_scf = True
+        else:
+            do_init_scf = self.init_scf
+            oqp.library.ints_1e(self.mol)
+
+        # compute energy
+        energies = self.sp.energy(do_init_scf=do_init_scf)
+
+        # compute QM/MM gradient
+        current_xyz = self.mol.get_system().reshape((-1, 3))
+        gradient_qm,gradient_mm=qmmm.openmm_gradient(current_xyz,self.mol.data["OQP::partial_charges"])
+        self.mol.data["OQP::mm_gradient"]=np.transpose(gradient_qm).tolist()
+
+        self.grad.grads = [self.istate]
+        grads = self.grad.gradient()
+        self.mol.energies = energies
+        self.mol.grads = grads
+
+        qmmm.gradient_qmmm=qmmm.form_gradient_qmmm(grads,gradient_mm)
+
+        # flatten data
+        energy = energies[self.istate]
+        grad = qmmm.gradient_qmmm.reshape(-1)
+
+        # evaluate metrics
+        de = energy - self.pre_energy
+        rmsd_step = np.mean((coordinates - self.pre_coord) ** 2) ** 0.5
+        max_step = np.amax(np.abs(coordinates - self.pre_coord))
+        rmsd_grad = np.mean(grad ** 2) ** 0.5
+        max_grad = np.amax(np.abs(grad))
+        self.metrics['itr'] = self.itr
+        self.metrics['de'] = de
+        self.metrics['rmsd_step'] = rmsd_step
+        self.metrics['max_step'] = max_step
+        self.metrics['rmsd_grad'] = rmsd_grad
+        self.metrics['max_grad'] = max_grad
+
+        # store energy and coordinates
+        self.pre_energy = energy
+        self.pre_coord = coordinates.copy()
+        dump_data((self.itr, self.atoms, coordinates, energy, de, rmsd_step, max_step, rmsd_grad, max_grad),
+                  title='OPTIMIZATION', fpath=self.mol.log_path)
+
+        return energy, grad.reshape((self.natom, 3))
+

pyoqp/oqp/library/libopenmm.py

@@ -0,0 +1,128 @@
+"""OQP OpenMM MD class"""
+import copy
+import oqp
+import numpy as np
+import scipy as sc
+from oqp.library.single_point import SinglePoint, Gradient, LastStep
+from oqp.utils.file_utils import dump_log, dump_data
+import oqp.utils.qmmm as qmmm 
+import openmm as mm
+import openmm.app as app
+import openmm.unit as unit
+
+
+class QMMM_MD:
+    def __init__(self, mol):
+        self.mol = mol
+        self.nSteps = mol.config['qmmm']['nsteps']
+        self.timeStep = mol.config['qmmm']['timestep']
+        self.init_scf = mol.config['optimize']['init_scf']
+        self.nstate = mol.config['tdhf']['nstate']
+        self.natom = mol.data['natom']
+        self.atoms = mol.get_atoms().reshape((self.natom, 1))
+        self.sp = SinglePoint(mol)
+        self.grad = Gradient(mol)
+
+        dump_log(self.mol, title='PyOQP: Entering QM/MM molecular dynamics')
+
+    def run_md(self):
+
+        pdb = qmmm.pdb0
+
+#Create an empty system and add particles based on PDB topology
+        system = mm.System()
+        for atom in pdb.topology.atoms():
+           system.addParticle(atom.element.mass)
+        number_of_particles = system.getNumParticles()
+
+#Add QM/MM gradient and energy to system as a CustomExternalForce
+        oqp_qmmm = mm.CustomExternalForce('grad_x*x + grad_y*y + grad_z*z + qmmm_energy - ecorr')
+        force_index = system.addForce(oqp_qmmm)
+        oqp_qmmm.addPerParticleParameter('grad_x')
+        oqp_qmmm.addPerParticleParameter('grad_y')
+        oqp_qmmm.addPerParticleParameter('grad_z')
+
+#First call to OQP and define initial parameters as QM/MM energy and gradient
+        qmmm_energy, qmmm_grad = self.oqp_qmmm_single_point(0)
+# Energy
+        qmmm_energy = (qmmm_energy/number_of_particles/qmmm.kj_per_mol_to_hartree)*unit.kilojoules_per_mole
+        oqp_qmmm.addGlobalParameter('qmmm_energy', qmmm_energy)
+# Gradient 
+        qmmm_grad = (qmmm_grad/qmmm.kj_per_mol_to_hartree/qmmm.bohr_to_nm) * unit.kilojoules_per_mole/unit.nanometer
+        ecorr=0.0*unit.kilojoules_per_mole
+        for i in range(number_of_particles):
+            oqp_qmmm.addParticle(i, qmmm_grad[i])
+            ecorr+=pdb.positions[i][0]*qmmm_grad[i][0]
+            ecorr+=pdb.positions[i][1]*qmmm_grad[i][1]
+            ecorr+=pdb.positions[i][2]*qmmm_grad[i][2]
+        ecorr/=number_of_particles
+#This is an energy correction to eliminate the grad_x*x + grad_y*y + grad_z*z contribution
+        oqp_qmmm.addGlobalParameter('ecorr', ecorr)
+
+#Create simulation integrator
+        integrator = mm.VerletIntegrator(self.timeStep*unit.femtoseconds)
+        simulation = app.Simulation(pdb.topology, system, integrator)
+        simulation.context.setPositions(pdb.positions)
+        simulation.context.setVelocitiesToTemperature(300.0*unit.kelvin)
+
+#Add reporters
+        simulation.reporters.append(app.PDBReporter('trajectory.pdb', 1))
+        oqp_output=open(self.mol.log+'.md','w')
+        oqp_output.write(f"\n\nStarting NVE dynamics using Verlet algorithm:\n")
+        simulation.reporters.append(app.StateDataReporter(oqp_output, 1, step=True, time=True, separator=" ",totalEnergy=True, kineticEnergy=True, potentialEnergy=True, temperature=False, speed=True,elapsedTime=True))
+
+        state = simulation.context.getState(getPositions=True,getForces=True)
+
+        for step in range(self.nSteps):
+
+            simulation.step(1)
+
+#Get the current positions from simulation context and update the pdb
+            state = simulation.context.getState(getPositions=True)
+            qmmm.pdb0.positions = state.getPositions()[:number_of_particles]
+
+#Perform OQP calculation and update parameters in context 
+            qmmm_energy, qmmm_grad = self.oqp_qmmm_single_point(step+1)
+            qmmm_energy = qmmm_energy/number_of_particles*unit.kilojoules_per_mole/qmmm.kj_per_mol_to_hartree
+            simulation.context.setParameter('qmmm_energy', qmmm_energy)
+
+            qmmm_grad = (qmmm_grad/qmmm.kj_per_mol_to_hartree/qmmm.bohr_to_nm) * unit.kilojoules_per_mole/unit.nanometer
+            ecorr=0.0*unit.kilojoules_per_mole
+            for i in range(number_of_particles):
+                oqp_qmmm.setParticleParameters(i, i, qmmm_grad[i])
+                ecorr+=pdb.positions[i][0]*qmmm_grad[i][0]
+                ecorr+=pdb.positions[i][1]*qmmm_grad[i][1]
+                ecorr+=pdb.positions[i][2]*qmmm_grad[i][2]
+            ecorr/=number_of_particles
+            simulation.context.setParameter('ecorr', ecorr)
+            oqp_qmmm.updateParametersInContext(simulation.context)
+
+    def oqp_qmmm_single_point(self,itr):
+
+        dump_log(self.mol, title='PyOQP: QM/MM molecular dynamics (NVE) using Verlet algorithm ')
+
+        num_atoms, x, y, z, q, mass = qmmm.openmm_system()
+        coordinates_qm=np.array(x + y + z).reshape((3, num_atoms)).T.reshape(-1)
+
+        self.mol.update_system(coordinates_qm)
+
+        # compute 1e integral
+        if itr == 0:
+            do_init_scf = True
+        else:
+            do_init_scf = self.init_scf
+            oqp.library.ints_1e(self.mol)
+
+        # compute energy
+        energies = self.sp.energy(do_init_scf=do_init_scf)
+
+        # compute QM/MM gradient
+        grads = self.grad.gradient()
+
+        # flatten data
+        energy = energies[qmmm.istate]
+        grad = qmmm.gradient_qmmm
+
+        return energy, grad
+
+