Reuse ForceField Parameters From Previous Calculations

ForceField Technical Workflow automation

Reuse force field atom types and charges assigned in a previous AMS calculation run by using the get_forcefield_params method.

Downloads: Notebook | Script ?

Initial Imports

import sys
from scm.plams import AMSJob, Settings, from_smiles

Run/Load Job with ForceField Information

First run a reference calculation where charges are guessed (using DFTB by default):

ref_job = AMSJob.from_input(
    """
Task GeometryOptimization

GeometryOptimization
   Convergence Step=1.0e-3
End

System
   Atoms
      C 0.0 0.0 0.0
      O 1.13 0.0 0.0
      C 0.0 0.0 2.1
      O 1.13 0.0 1.9
   End
End

Engine ForceField
   Verbosity Verbose
   GuessCharges True
EndEngine
"""
)

# Alternatively, load a previously run calculation
# ref_job = AMSJob.load_external("./plams_workdir/plamsjob/ams.rkf")

ref_job.run();

[25.03|09:47:58] JOB plamsjob STARTED
[25.03|09:47:58] JOB plamsjob RUNNING
[25.03|09:47:59] JOB plamsjob FINISHED
[25.03|09:47:59] JOB plamsjob SUCCESSFUL

Reuse ForceField Parameters

Extract the charges and types from the job results and add them as properties on the molecule:

charges, types, patch = ref_job.results.get_forcefield_params()

mol = ref_job.molecule[""].copy()

for i, at in enumerate(mol.atoms):
    at.properties.ForceField.Charge = charges[i]
    at.properties.ForceField.Type = types[i]

sett = Settings()
sett.input.AMS.Task = "SinglePoint"
sett.input.ForceField.Type = "UFF"

# Create a patch file if required
if patch:
    with open("patch.dat", "w") as outfile:
        outfile.write(str(patch))
        outfile.close()
    # For example with:
    # sett.input.ForceField.GAFF.ForceFieldPatchFile = "patch.dat"

job = AMSJob(molecule=mol, settings=sett)

print(job.get_input())

Task SinglePoint

System
  Atoms
              C       0.0000000000       0.0000000000       0.0000000000 ForceField.Charge=0.2881959744167275 ForceField.Type=C_1
              O       1.1300000000       0.0000000000       0.0000000000 ForceField.Charge=-0.2676126103828702 ForceField.Type=O_2
              C       0.0000000000       0.0000000000       2.1000000000 ForceField.Charge=0.2536150412119178 ForceField.Type=C_1
              O       1.1300000000       0.0000000000       1.9000000000 ForceField.Charge=-0.27419840524497996 ForceField.Type=O_2
  End
End

Engine ForceField
  Type UFF
EndEngine

job.run();

[25.03|09:47:59] JOB plamsjob STARTED
[25.03|09:47:59] Renaming job plamsjob to plamsjob.002
[25.03|09:47:59] JOB plamsjob.002 RUNNING
[25.03|09:48:00] JOB plamsjob.002 FINISHED
[25.03|09:48:00] JOB plamsjob.002 SUCCESSFUL

Python Script

#!/usr/bin/env python
# coding: utf-8

# ## Initial Imports

import sys
from scm.plams import AMSJob, Settings, from_smiles


# ## Run/Load Job with ForceField Information

# First run a reference calculation where charges are guessed (using DFTB by default):

ref_job = AMSJob.from_input(
    """
Task GeometryOptimization

GeometryOptimization
   Convergence Step=1.0e-3
End

System
   Atoms
      C 0.0 0.0 0.0
      O 1.13 0.0 0.0
      C 0.0 0.0 2.1
      O 1.13 0.0 1.9
   End
End

Engine ForceField
   Verbosity Verbose
   GuessCharges True
EndEngine
"""
)

# Alternatively, load a previously run calculation
# ref_job = AMSJob.load_external("./plams_workdir/plamsjob/ams.rkf")

ref_job.run()
# ## Reuse ForceField Parameters

# Extract the charges and types from the job results and add them as properties on the molecule:

charges, types, patch = ref_job.results.get_forcefield_params()


mol = ref_job.molecule[""].copy()

for i, at in enumerate(mol.atoms):
    at.properties.ForceField.Charge = charges[i]
    at.properties.ForceField.Type = types[i]


sett = Settings()
sett.input.AMS.Task = "SinglePoint"
sett.input.ForceField.Type = "UFF"


# Create a patch file if required
if patch:
    with open("patch.dat", "w") as outfile:
        outfile.write(str(patch))
        outfile.close()
    # For example with:
    # sett.input.ForceField.GAFF.ForceFieldPatchFile = "patch.dat"


job = AMSJob(molecule=mol, settings=sett)


print(job.get_input())


job.run()