Quickstart#
First let import the modules to work with. Some functionalities of TOFF may not be available on PyPi. This could be because the new version is not yet release, but in brief it is going to be. You could install directly from the repo if some problem pops up. Just paste the following in a code cell:
! pip install git+https://github.com/ale94mleon/TOFF.git@main
As module#
from toff import Parameterize
from rdkit import Chem
import tempfile, os, yaml
# We will use this temporal directory, you could change it in your case to, for example wd = '.'
tmp_dir = tempfile.TemporaryDirectory()
wd = tmp_dir.name
Initialize the class#
The first steep is to create an instance of the class Parameterize.init (follow the link to know all possible arguments to initialize the class). For now we will use the default parameters, we will only specify the out_dir.
parameterizer = Parameterize(out_dir=wd)
parameterizer
Parameterize(force_field_code = openff_unconstrained-2.0.0.offxml, ext_types = [top pdb gro], hmr_factor = None, overwrite = False, safe_naming_prefix = None, out_dir = /tmp/tmpsn2h1yw_)
call the class#
You could see what ara all the arguments to specify here Parameterize.call.
parameterizer(Chem.MolFromSmiles('CCO'))
os.listdir(wd)
['MOL.gro', 'MOL.pdb', 'MOL.top']
Take a look on the .top file#
with open(os.path.join(wd, 'MOL.top'), 'r') as f:
text = f.read()
print(text)
;
; File /tmp/tmpsn2h1yw_/MOL.top was generated
; By user: ale (1000)
; On host:ale-MS-7B86
; At date:Tue. May 2 11:49:36 2023
;
; This is a standalone topology file
;
; Created by:
; ParmEd: ipykernel_launcher.py, VERSION4.1.0
; Executable: ipykernel_launcher.py
; Library dir: /usr/local/gromacs/share/gromacs/top
; Command line:
; /home/ale/anaconda3/envs/toff/lib/python3.10/site-packages/ipykernel_launcher.py --ip=127.0.0.1 --stdin=9003 --control=9001 --hb=9000 --Session.signature_scheme="hmac-sha256" --Session.key=b"e47eeab4-8252-4f27-97b3-94bb22f4ec2b" --shell=9002 --transport="tcp" --iopub=9004 --f=/home/ale/.local/share/jupyter/runtime/kernel-v2-112377m8T2Gm0VzHsA.json
;
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 no 1 0.83333333
[ atomtypes ]
; name at.num mass charge ptype sigma epsilon
C1 6 12.010780 0.00000000 A 0.33795318 0.45538912
O1 8 15.999430 0.00000000 A 0.299716 0.87643726
H1 1 1.007947 0.00000000 A 0.26445434 0.066021356
H2 1 1.007947 0.00000000 A 0.25832257 0.068656285
H3 1 1.007947 0.00000000 A 0.053453923 5.1571983e-05
[ moleculetype ]
; Name nrexcl
MOL 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 MOL rtp MOL q 0.0
1 C1 1 MOL C1 1 -0.13610010 12.010780 ; qtot -0.136100
2 C1 1 MOL C2 2 0.12639990 12.010780 ; qtot -0.009700
3 O1 1 MOL O1 3 -0.59980010 15.999430 ; qtot -0.609500
4 H1 1 MOL H1 4 0.04236690 1.007947 ; qtot -0.567133
5 H1 1 MOL H2 5 0.04236690 1.007947 ; qtot -0.524766
6 H1 1 MOL H3 6 0.04236690 1.007947 ; qtot -0.482400
7 H2 1 MOL H4 7 0.04319990 1.007947 ; qtot -0.439200
8 H2 1 MOL H5 8 0.04319990 1.007947 ; qtot -0.396000
9 H3 1 MOL H6 9 0.39599990 1.007947 ; qtot 0.000000
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1 0.15219 221435.259290
1 4 1 0.10939 309655.084322
1 5 1 0.10939 309655.084322
1 6 1 0.10939 309655.084322
2 3 1 0.14273 276118.879749
2 7 1 0.10939 309655.084322
2 8 1 0.10939 309655.084322
3 9 1 0.09717 454823.212172
[ pairs ]
; ai aj funct c0 c1 c2 c3
3 4 1 0.282085170 0.120274245
3 5 1 0.282085170 0.120274245
3 6 1 0.282085170 0.120274245
4 7 1 0.261388456 0.033662966
5 7 1 0.261388456 0.033662966
6 7 1 0.261388456 0.033662966
4 8 1 0.261388456 0.033662966
5 8 1 0.261388456 0.033662966
6 8 1 0.261388456 0.033662966
1 9 1 0.195703550 0.002423083
7 9 1 0.155888247 0.000940843
8 9 1 0.155888247 0.000940843
[ angles ]
; ai aj ak funct c0 c1 c2 c3
1 2 3 1 116.5475863 445.222087
1 2 7 1 116.5475863 445.222087
1 2 8 1 116.5475863 445.222087
2 1 4 1 116.5475863 445.222087
2 1 5 1 116.5475863 445.222087
2 1 6 1 116.5475863 445.222087
2 3 9 1 110.3538806 544.678275
3 2 7 1 116.5475863 445.222087
3 2 8 1 116.5475863 445.222087
4 1 5 1 115.6031000 408.161690
4 1 6 1 115.6031000 408.161690
5 1 6 1 115.6031000 408.161690
7 2 8 1 115.6031000 408.161690
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
1 2 3 9 1 0.0000000 1.4415479 3
1 2 3 9 1 0.0000000 0.5709304 1
3 2 1 4 1 0.0000000 0.4667359 3
3 2 1 4 1 0.0000000 1.4361539 1
3 2 1 5 1 0.0000000 0.4667359 3
3 2 1 5 1 0.0000000 1.4361539 1
3 2 1 6 1 0.0000000 0.4667359 3
3 2 1 6 1 0.0000000 1.4361539 1
4 1 2 7 1 0.0000000 0.7999501 3
4 1 2 8 1 0.0000000 0.7999501 3
5 1 2 7 1 0.0000000 0.7999501 3
5 1 2 8 1 0.0000000 0.7999501 3
6 1 2 7 1 0.0000000 0.7999501 3
6 1 2 8 1 0.0000000 0.7999501 3
7 2 3 9 1 0.0000000 1.2662416 3
8 2 3 9 1 0.0000000 1.2662416 3
[ system ]
; Name
Generic title
[ molecules ]
; Compound #mols
MOL 1
We could reuse the instance of Parameterize and used in other molecule
parameterizer(
input_mol=Chem.MolFromSmiles('NCCCC'),
mol_resi_name='LIG',
)
os.listdir(wd)
Charges will be corrected: partial_charge - formal_charge = -1e-06
After correction: partial_charge - formal_charge = 0.0.
['LIG.top', 'MOL.gro', 'MOL.pdb', 'LIG.pdb', 'MOL.top', 'LIG.gro']
We would like to use Hydrogen Mass Repartitioning. This allows to increases the integration time step from 2 fs (usually) to 4fs. Here we have to create a new instance of Parameterize and set the hmr_factor. We will use 3, which means that the mass of the hydrogens atoms will increase in a factor of 3 using the mass of the heavy atoms attached to them. We also specify a safe_naming_prefix, this avoid possible naming conflicts with pother force fields.
parameterizer_hmr = Parameterize(out_dir=wd, hmr_factor=3, safe_naming_prefix='z')
parameterizer_hmr
Parameterize(force_field_code = openff_unconstrained-2.0.0.offxml, ext_types = [top pdb gro], hmr_factor = 3, overwrite = False, safe_naming_prefix = z, out_dir = /tmp/tmpsn2h1yw_)
parameterizer_hmr(
input_mol=Chem.MolFromSmiles('CCF'),
mol_resi_name='HMR',
)
os.listdir(wd)
['LIG.top',
'MOL.gro',
'HMR.top',
'HMR.gro',
'MOL.pdb',
'LIG.pdb',
'MOL.top',
'LIG.gro',
'HMR.pdb']
with open(os.path.join(wd, 'HMR.top'), 'r') as f:
text = f.read()
print(text)
;
; File /tmp/tmpsn2h1yw_/HMR.top was generated
; By user: ale (1000)
; On host:ale-MS-7B86
; At date:Tue. May 2 11:49:37 2023
;
; This is a standalone topology file
;
; Created by:
; ParmEd: ipykernel_launcher.py, VERSION4.1.0
; Executable: ipykernel_launcher.py
; Library dir: /usr/local/gromacs/share/gromacs/top
; Command line:
; /home/ale/anaconda3/envs/toff/lib/python3.10/site-packages/ipykernel_launcher.py --ip=127.0.0.1 --stdin=9003 --control=9001 --hb=9000 --Session.signature_scheme="hmac-sha256" --Session.key=b"e47eeab4-8252-4f27-97b3-94bb22f4ec2b" --shell=9002 --transport="tcp" --iopub=9004 --f=/home/ale/.local/share/jupyter/runtime/kernel-v2-112377m8T2Gm0VzHsA.json
;
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 no 1 0.83333333
[ atomtypes ]
; name at.num mass charge ptype sigma epsilon
zC1 6 12.010780 0.00000000 A 0.33795318 0.45538912
zF1 9 18.998403 0.00000000 A 0.31181455 0.255224
zH1 1 1.007947 0.00000000 A 0.26445434 0.066021356
zH2 1 1.007947 0.00000000 A 0.25832257 0.068656285
[ moleculetype ]
; Name nrexcl
HMR 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 HMR rtp HMR q 0.0
1 zC1 1 HMR C1 1 -0.13097500 6.034621 ; qtot -0.130975
2 zC1 1 HMR C2 2 0.15902500 8.026674 ; qtot 0.028050
3 zF1 1 HMR F1 3 -0.24917500 18.998403 ; qtot -0.221125
4 zH1 1 HMR H1 4 0.04782500 3.000000 ; qtot -0.173300
5 zH1 1 HMR H2 5 0.04782500 3.000000 ; qtot -0.125475
6 zH1 1 HMR H3 6 0.04782500 3.000000 ; qtot -0.077650
7 zH2 1 HMR H4 7 0.03882500 3.000000 ; qtot -0.038825
8 zH2 1 HMR H5 8 0.03882500 3.000000 ; qtot 0.000000
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1 0.15219 221435.259290
1 4 1 0.10939 309655.084322
1 5 1 0.10939 309655.084322
1 6 1 0.10939 309655.084322
2 3 1 0.13587 284889.572611
2 7 1 0.10939 309655.084322
2 8 1 0.10939 309655.084322
[ pairs ]
; ai aj funct c0 c1 c2 c3
3 4 1 0.288134446 0.064904227
3 5 1 0.288134446 0.064904227
3 6 1 0.288134446 0.064904227
4 7 1 0.261388456 0.033662966
5 7 1 0.261388456 0.033662966
6 7 1 0.261388456 0.033662966
4 8 1 0.261388456 0.033662966
5 8 1 0.261388456 0.033662966
6 8 1 0.261388456 0.033662966
[ angles ]
; ai aj ak funct c0 c1 c2 c3
1 2 3 1 116.5475863 445.222087
1 2 7 1 116.5475863 445.222087
1 2 8 1 116.5475863 445.222087
2 1 4 1 116.5475863 445.222087
2 1 5 1 116.5475863 445.222087
2 1 6 1 116.5475863 445.222087
3 2 7 1 116.5475863 445.222087
3 2 8 1 116.5475863 445.222087
4 1 5 1 115.6031000 408.161690
4 1 6 1 115.6031000 408.161690
5 1 6 1 115.6031000 408.161690
7 2 8 1 115.6031000 408.161690
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
3 2 1 4 1 0.0000000 0.4189792 3
3 2 1 4 1 0.0000000 1.4026079 1
3 2 1 5 1 0.0000000 0.4189792 3
3 2 1 5 1 0.0000000 1.4026079 1
3 2 1 6 1 0.0000000 0.4189792 3
3 2 1 6 1 0.0000000 1.4026079 1
4 1 2 7 1 0.0000000 0.7999501 3
4 1 2 8 1 0.0000000 0.7999501 3
5 1 2 7 1 0.0000000 0.7999501 3
5 1 2 8 1 0.0000000 0.7999501 3
6 1 2 7 1 0.0000000 0.7999501 3
6 1 2 8 1 0.0000000 0.7999501 3
[ system ]
; Name
Generic title
[ molecules ]
; Compound #mols
HMR 1
As command line#
We need first to create a yaml configuration file that is nothing else that a more human-readable python dictionary. The content of this file is all the variables needed for initialize and call the parameterize class. The only mandatory argument is input_mol and in this case must be a path to a file with the following extensions: smi, inchi, mol, mol2. In the case of smi and inchi only will be redded the first line of the file and must be a valid SMILES or InChi strings respectively.
# First let save a molecule
Chem.MolToMolFile(Chem.MolFromSmiles('CCCBr'), os.path.join(wd, 'my_file.mol'))
os.listdir(wd)
['LIG.top',
'MOL.gro',
'HMR.top',
'HMR.gro',
'MOL.pdb',
'LIG.pdb',
'MOL.top',
'my_file.mol',
'LIG.gro',
'HMR.pdb']
config = {
'force_field_code': 'openff-1.3.1.offxml', # by default is openff_unconstrained-2.0.0.offxml
'ext_types': ['mol3', 'top'], # by default is ['top', 'pdb', 'gro']. And ypou could pick from: ‘pdb’, ‘pqr’, ‘cif’,’pdbx’, ‘parm7’, ‘prmtop’, ‘psf’, ‘top’, ‘gro’, ‘mol2’, ‘mol3’, ‘crd’, ‘rst7’, ‘inpcrd’, ‘restrt’, ‘ncrst’
'hmr_factor': 2.5, # by default None
'overwrite': True, # by default False
'out_dir': os.path.join(wd, 'cmd'), # by default '.'
'input_mol': os.path.join(wd, 'my_file.mol'), # The only parameter that MUST BE PROVIDED!
'mol_resi_name': 'CMD', # by default 'MOL'
}
with open(os.path.join(wd, 'config.yml'), 'w') as cf:
yaml.dump(config, cf)
print(f"This is how it looks the wd: {os.listdir(wd)}")
with open(os.path.join(wd, 'config.yml'), 'r') as cfy:
text = cfy.read()
print("===========config.yml=======================")
print(text)
print("============================================")
# let's also get the path to the config.yml file to passed to the command line.
print(f"Path to the config file: {os.path.join(wd, 'config.yml')}")
This is how it looks the wd: ['LIG.top', 'MOL.gro', 'HMR.top', 'HMR.gro', 'MOL.pdb', 'LIG.pdb', 'MOL.top', 'my_file.mol', 'LIG.gro', 'config.yml', 'HMR.pdb']
===========config.yml=======================
ext_types:
- mol3
- top
force_field_code: openff-1.3.1.offxml
hmr_factor: 2.5
input_mol: /tmp/tmpsn2h1yw_/my_file.mol
mol_resi_name: CMD
out_dir: /tmp/tmpsn2h1yw_/cmd
overwrite: true
============================================
Path to the config file: /tmp/tmpsn2h1yw_/config.yml
# Change the path properly. Yur temporal path is printed in the cell before and for sure it will not be the same as mine.
! parameterize /tmp/tmpsn2h1yw_/config.yml
You are using toff:0.0.2
Parameterize(force_field_code = openff-1.3.1.offxml, ext_types = [mol3 top], hmr_factor = 2.5, overwrite = True, safe_naming_prefix = None, out_dir = /tmp/tmpsn2h1yw_/cmd)
# Check if the file were saved
os.listdir(os.path.join(wd, 'cmd'))
['CMD.top', 'CMD.mol3']