Example: Loading MM charges for regions
***************************************

In this example we consider an OH- with an H3O+ fragment. As the
charges on the fragments are kept fixed, the formation of two water
molecules is avoided.

First we "estimate" the charges for the two fragments with a DFTB
calculation.

These charges are then loaded for the correct regions in the total
system. Observe that this is done in the "System" block, see the
System definition section of the AMS manual.

We do this first for a QUILD-like setup (mechanical embedding), and
next for a QMMM calculation with electrostatic coupling.

Note that this example makes use of the python scripting tools PLAMS
and PISA, available in the AMS Python Stack.

"Download LoadCharges.py"

   #!/usr/bin/env amspython

   from scm.base import ChemicalSystem
   from scm.input_classes import drivers, engines
   from scm.plams import AMSJob, JobRunner, config, Settings
   import copy

   # Here we treat H3O+ as qm and OH- as the MM region (Optimizing without regions gives two H2O molecules)
   # We do this with a QUILD setup (mechanical embedding) and electrostatic embedding (QMMM)
   # We obtain the charges from a DFTB calculation
   # In this case the results (QUILD vs. QMMM) are quite similar as apparently the OH does not polarize the QM region much

   # First we do two DFTB calculations on the two fragments
   config.default_jobrunner = JobRunner(parallel=True)
   settings = Settings()
   settings.runscript.nproc = 1

   h30 = ChemicalSystem(
       """
   System
       Atoms
           O -1.527946410885647 -0.2107366711137158 -0.0008116899510243671
           H -0.8459142126057956 0.3517312394359257 0.4094504676540848
           H -1.834953147575289 0.1051014241823828 -0.8704652381864062
           H -1.328032016244278 -1.164422847242489 0.02894848344144469
       End
       Charge 1.0
       GuessBonds True
   End
   """
   )

   oh = ChemicalSystem(
       """
   System
       Atoms
           O 0.6370858511871781 -0.3378071707560572 -0.0006181020627287671
           H 1.318474396634582 0.2241299231185073 0.4092568796869673
       End
       Charge -1.0
       GuessBonds True
   End
   """
   )

   fragments = {"H3O+": h30, "OH-": oh}

   results_dftb = {}
   for name, mol in fragments.items():
       driver_dftb = drivers.AMS()
       driver_dftb.Task = "SinglePoint"
       driver_dftb.Properties.Charges = True
       driver_dftb.GeometryOptimization.Convergence.Gradients = 1.0e-6
       driver_dftb.Engine = engines.DFTB()
       job = AMSJob(molecule=mol, settings=settings + Settings({"input": driver_dftb}), name=name)
       results_dftb[name] = job.run()

   # Now we run it in a QUILD-like setup (mechanical embedding) and in a QMMM-like setup
   h20_2 = ChemicalSystem(
       """
   System
       Atoms
           O -1.527946410885647 -0.2107366711137158 -0.0008116899510243671 region=QM
           H -0.8459142126057956 0.3517312394359257 0.4094504676540848     region=QM
           H -1.834953147575289 0.1051014241823828 -0.8704652381864062     region=QM
           H -1.328032016244278 -1.164422847242489 0.02894848344144469     region=QM
           O 0.6370858511871781 -0.3378071707560572 -0.0006181020627287671 region=MM
           H 1.318474396634582 0.2241299231185073 0.4092568796869673       region=MM
       End
       GuessBonds True
   End
   """
   )

   qm_charges = results_dftb["H3O+"].get_charges()
   mm_charges = results_dftb["OH-"].get_charges()
   h20_2.enable_atom_attributes("forcefield")

   for i, charge in zip(h20_2.get_atoms_in_region("QM"), qm_charges):
       h20_2.atoms[i].forcefield.charge = charge

   for i, charge in zip(h20_2.get_atoms_in_region("MM"), mm_charges):
       h20_2.atoms[i].forcefield.charge = charge

   results = {}
   for name in ["quild", "qmmm"]:
       driver = copy.deepcopy(driver_dftb)
       driver.Task = "GeometryOptimization"
       driver.Engine = engines.Hybrid()

       driver.Engine.Engine[0] = engines.DFTB()
       driver.Engine.Engine[0].Model = "GFN1-xTB"
       driver.Engine.Engine[1] = engines.ForceField()

       if name == "quild":
           driver.Engine.Energy.Term[0].Factor = 1.0
           driver.Engine.Energy.Term[0].Region = "*"
           driver.Engine.Energy.Term[0].EngineID = "ForceField"

           driver.Engine.Energy.Term[1].Factor = -1.0
           driver.Engine.Energy.Term[1].Region = "QM"
           driver.Engine.Energy.Term[1].EngineID = "ForceField"
           driver.Engine.Energy.Term[1].Charge = 1.0

           driver.Engine.Energy.Term[2].Factor = 1.0
           driver.Engine.Energy.Term[2].Region = "QM"
           driver.Engine.Energy.Term[2].EngineID = "DFTB"
           driver.Engine.Energy.Term[2].Charge = 1.0

       elif name == "qmmm":
           driver.Engine.QMMM.QMRegion = "QM"
           driver.Engine.QMMM.QMEngineID = "DFTB"
           driver.Engine.QMMM.MMEngineID = "ForceField"
           driver.Engine.QMMM.QMCharge = 1.0
           driver.Engine.QMMM.MMCharge = -1.0

       job = AMSJob(molecule=h20_2, settings=settings + Settings({"input": driver}), name=name)
       results[name] = job.run()

   results = {k: (v.get_charges(), v.get_main_system()) for k, v in results.items()}

   print("start of report")
   print("method distance charges")
   for name, (charges, system) in results.items():
       print(
           f"{name} {system.get_distance(0, 4):5.3f} {charges[0]:5.3f} {charges[1]:5.3f} {charges[2]:5.3f} {charges[3]:5.3f} {charges[4]:5.3f} {charges[5]:5.3f}"
       )
   print("end of report")
