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.

Download LoadCharges.run

#!/bin/sh

# 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

report=report.txt

echo "method distance charges" > $report


# first we do two DFTB calculations on the two fragments

export AMS_JOBNAME=H2O+.dftb

rm -rf $AMS_JOBNAME.results

"$AMSBIN/ams" << eor

Task SinglePoint

Properties Charges=yes

GeometryOptimization
  Convergence Gradients=1.0e-6
End

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

Engine DFTB
EndEngine
    

eor

export AMS_JOBNAME=OH-.dftb

rm -rf $AMS_JOBNAME.results

"$AMSBIN/ams" << eor

Task SinglePoint

Properties Charges=yes

GeometryOptimization
  Convergence Gradients=1.0e-6
End

System
    Atoms
        O 0.6370858511871781 -0.3378071707560572 -0.0006181020627287671
        H 1.318474396634582 0.2241299231185073 0.4092568796869673      
    End
    Charge -1.0
    GuessBonds True
End

Engine DFTB
EndEngine
    

eor

# Now we run it in a QUILD-like setup (mechanical embedding)

export AMS_JOBNAME=quild

rm -rf $AMS_JOBNAME.results

"$AMSBIN/ams" << eor

Task GeometryOptimization

Properties Charges=yes

GeometryOptimization
  Convergence Gradients=1.0e-6
End

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

    LoadForceFieldCharges region=QM file=H2O+.dftb.results
    LoadForceFieldCharges region=MM file=OH-.dftb.results
End

Engine Hybrid
    Energy
       Term Factor=1.0  Region=*  EngineID=ForceField
       Term Factor=-1.0 Region=QM EngineID=ForceField     Charge=1.0
       Term Factor=1.0  Region=QM EngineID=DFTB           Charge=1.0
    End

    Engine DFTB
        Model GFN1-xTB
    EndEngine
    
    Engine ForceField
    EndEngine
EndEngine

eor


ddd=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#5`
eee=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k AMSResults%Charges#5.3f`

echo "quild  $charge $ddd $eee" >> $report

# Now we run it in a QMMM-like setup


export AMS_JOBNAME=qmmm

rm -rf $AMS_JOBNAME.results

"$AMSBIN/ams" << eor

Properties Charges=yes

Task GeometryOptimization

GeometryOptimization
  Convergence Gradients=1.0e-6
End

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

    LoadForceFieldCharges region=QM file=H2O+.dftb.results
    LoadForceFieldCharges region=MM file=OH-.dftb.results
End

Engine Hybrid
    QMMM QMRegion=QM QMEngineID=DFTB MMEngineID=ForceField QMCharge=1.0 MMCharge=-1.0

    Engine DFTB
        Model GFN1-xTB
    EndEngine
    
    Engine ForceField
    EndEngine
EndEngine

eor

ddd=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#5`
eee=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k AMSResults%Charges#5.3f`

echo "qmmm   $charge $ddd $eee" >> $report


echo "start of report"
cat $report
echo "end of report"