JMorado/ParaMol

Adaptive parameterization question

mmagithub opened this issue · 1 comments

Hi,
I am trying to optimize an existing ligand amber parameter topology file using the adaptive parameterization scheme (DFTB+). At the updated parameters writing step, I am observing that the saved updated amber prmtop file stays the same as the original prmtop file, even though the 'lig_adaptive_param.ff' got updated, i.e. it is slightly different from the file saved before the adaptive parameterization. I am wondering if I am missing something. Here is my script:

#########################
import time
from time import process_time
t1_start = process_time()

print(t1_start)

from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *
from ParaMol.QM_engines.qm_engine import *
from ParaMol.Tasks.adaptive_parametrization import *
from ParaMol.Utils.settings import *
from ParaMol.Utils.amber_symmetrizer import *

---------------------------------------------------------

Preparation

---------------------------------------------------------

Create the OpenMM engine for lig

openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='lig.prmtop', crd_file='lig.inpcrd')

Create ParaMol System

lig = ParaMolSystem(name="lig", engine=openmm_system, n_atoms=67, n_cpus=4)

Create ParaMol's force field representation and ask to parametrize bonds, angles and torsions

lig.force_field.create_force_field(opt_bonds=True, opt_angles=True, opt_torsions=True)

Create ParaMol settings instance

paramol_settings = Settings()

The objective function will contain a energy, force and regularization term

paramol_settings.properties["include_energies"] = True
paramol_settings.properties["include_forces"] = True
paramol_settings.properties["include_regularization"] = True

Create the ASE calculator

from ase.calculators.dftb import *

calc = Dftb(Hamiltonian_='DFTB', # line is included by default
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_H='s',
Hamiltonian_MaxAngularMomentum_O='p',
Hamiltonian_MaxAngularMomentum_C='p',
Hamiltonian_MaxAngularMomentum_N="p",
Hamiltonian_MaxAngularMomentum_S="p",
Hamiltonian_MaxAngularMomentum_Br='p',
Hamiltonian_MaxAngularMomentum_F='p',
Hamiltonian_MaxAngularMomentum_Cl='p',
Hamiltonian_MaxAngularMomentum_I='p',
Hamiltonian_MaxAngularMomentum_P='p',
Hamiltonian_Dispersion="DftD3 { \n s6=1.000 \n s8=0.5883 \n Damping = BeckeJohnson { \n a1=0.5719 \n a2=3.6017 \n } \n }",
Hamiltonian_SCC='Yes',
Hamiltonian_SCCTolerance=1e-8,slako_dir='~/dftbplus-20.2.1/slako/3ob/3ob-3-1/')

Alternative, we could set the calculator in the settings

paramol_settings.qm_engine["qm_engine"] = "ase"
paramol_settings.qm_engine["ase"]["calculator"] = calc

amber_symmetrizer = AmberSymmetrizer("lig.prmtop")
amber_symmetrizer.get_amber_symmetries(lig.force_field)

Write symmetrized Force-Field to file

lig.force_field.write_ff_file("lig_sym.ff")

---------------------------------------------------------

Adaptive Parametrization

---------------------------------------------------------

adaptive_parametrization = AdaptiveParametrization()
systems, parameter_space, objective_function, optimizer = adaptive_parametrization.run_task(paramol_settings, [lig], rmsd_tol=0.01, max_iter=50, structures_per_iter=100, )

Write ParaMol Force Field file with final parameters

lig.force_field.write_ff_file("lig_adaptive_param.ff")

Update AMBER symmetrizer with new parameters

amber_symmetrizer.update_term_types_parameters(parameter_space.optimizable_parameters)

Write AMBER topology file (.prmtop)

amber_symmetrizer.save_prmtop("lig_opt.prmtop")
amber_symmetrizer.save_frcmod("lig_opt.frcmod")

t1_stop = process_time()
print("Elapsed time:", t1_stop, t1_start)
print("Elapsed time during the whole program in hours:",(t1_stop-t1_start)/3600)
#########################

Any clue,

Thanks,
Marawan

Hi,

I think you are still using an old version of the code (probably 1.0.2?). Please update your code through conda or by pointing to this repo as I have mentioned in a previous Issue.

The save_prmtop does not exist in the current version (you'd use simply save) and the way the writing and symmetrization of the FF are done has changed considerably since then.

Due to time constraints, I do not give support to deprecated versions of ParaMol. Furthermore, be aware that you may run into unexpected behavior as that version was released during a time wherein the code was still undergoing intense development.

Best regards,
João