Tutorials¶
Tutorial 1 - Benzene¶
This tutorial will describe in detail the use of QuickFF by means of the command line scripts, to generate a force field for benzene. This tutorial assumes that the following input files are available (examples of such file are provided with the source code in the directory share/systems/benzene):
gaussian.fchkA Gaussian Formatted Checkpoint file containing the ab initio equilibrium geometry and ab initio Hessian in equilibrium. This file can be generated by performing a freq job in Gaussian using your level of theory and basis set of choice. The following line should be added prior to the route section (i.e. the lines starting with #):
%chk=gaussian.chk
When the calculation terminated succesfully, make the fchk file using the following command:
formchk gaussian.chk gaussian.fchk
gaussian_mbis.h5A HDF5 file containing the atomic charges. Such a file can, for example, be generated using HORTON. If you have HORTON 2.x installed, the following command will derive the atomic charges in this example from the wavefunction in the
gaussian.fchkfile:horton-wpart.py gaussian.fchk gaussian_mbis.h5 mbis --grid=ultrafine
To determine the path in the HDF5 file to these charges, one can use the h5dump command (part of the Linux package hdf5-tools):
h5dump -n gaussian_mbis.h5
The output of this command for the gaussian_mbis.h5 file provided in the directory
share/tuturials/benzeneis:HDF5 "gaussian_mbis.h5" { FILE_CONTENTS { group / dataset /cartesian_multipoles dataset /change dataset /charges dataset /core_charges dataset /history_charges dataset /history_propars dataset /niter dataset /populations dataset /propars dataset /pseudo_populations dataset /pure_multipoles dataset /radial_moments dataset /valence_charges dataset /valence_widths } }From this output, we can conclude that, for this example, MBIS charges can be found in the path
/charges.
Covalent force field without non-bonding terms¶
As a first example, we will derive a force field containing only covalent terms.
No electrostatic nor van der Waals interactions will be included in the force
field. Furthermore, we use the built-in feature to automatically determine atom types according to
the level low (for benzene it does not make any difference which level is
chosen). This can very easily be done using the qff.py script (documentation on the available options can be
accessed by qff.py --help):
qff.py --ffatypes=low --suffix=_noei gaussian.fchk
The script will dump all relevant information to the screen, for this tutorial, the output is as follows:
________________/\\\_________/\\\\\\\\\\\\\\\__/\\\\\\\\\\\\\\\_________________
______________/\\\\/\\\\_____\/\\\///////////__\/\\\///////////_________________
_____________/\\\//\////\\\___\/\\\_____________\/\\\___________________________
_____________/\\\______\//\\\__\/\\\\\\\\\\\_____\/\\\\\\\\\\\__________________
_____________\//\\\______/\\\___\/\\\///////______\/\\\///////__________________
_______________\///\\\\/\\\\/____\/\\\_____________\/\\\________________________
__________________\////\\\//______\/\\\_____________\/\\\_______________________
______________________\///\\\\\\___\/\\\_____________\/\\\______________________
_________________________\//////____\///______________\///______________________
Welcom to QuickFF
a Python package to quickly derive force fields from ab initio input data
Written by
Louis Vanduyfhuys(1)*, Steven Vandenbrande(1) and Toon Verstraelen(1)
(1) Center for Molecular Modeling, Ghent University Belgium.
* mailto: Louis.Vanduyfhuys@UGent.be
USER louis
MACHINE Linux T085 4.15.0-39-generic #42-Ubuntu SMP Tue Oct 23
MACHINE 15:48:01 UTC 2018 x86_64
TIME 2018-11-29 16:23:09.283297
QUICKFF VERSION 2.2.1
PYTHON VERSION 3.6.3 |Anaconda custom (64-bit)| (default, Oct 13 2017,
PYTHON VERSION 12:02:49) [GCC 7.2.0]
NUMPY VERSION 1.13.3
SCIPY VERSION 0.19.1
MATPLOTLIB VERSION 2.1.0
CURRENT DIR /home/louis/build/quickff/share/tutorials/benzene
COMMAND LINE /home/louis/build/anaconda3/bin/qff.py --ffatypes=low
COMMAND LINE --suffix=_noei gaussian.fchk
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
INIT Initializing system
INIT Initializing Second order taylor reference for ai
INIT Initializing program
VAL Initializing valence force field
VAL Added 12 bond terms
VAL Added 18 bend terms (an)harmonic in the angle, 0 bend terms with
VAL 1+cos(angle) potential and 0 bend terms with 1-cos(4*theta) potential.
VAL Added 6 Harmonic and 0 SquareHarmonic out-of-plane distance terms
EQSET Resetting system coordinates to ab initio ref
EQSET Setting rest values to AI equilibrium values for tasks EQ_RV
EQSET Averaging force field parameters over master and slaves
PTGEN Resetting system coordinates to ab initio ref
PTGEN Constructing trajectories
PTEST Resetting system coordinates to ab initio ref
PTEST Estimating FF parameters from perturbation trajectories
PTPOST Averaging force field parameters over master and slaves
VAL Resetting system coordinates to ab initio ref
HCEST Resetting system coordinates to ab initio ref
HCEST Estimating force constants from Hessian cost for tasks HC_FC_DIAG HC_FC_
HCEST CROSS_ASS HC_FC_CROSS_ASA
PTEST Resetting system coordinates to ab initio ref
PTEST Estimating FF parameters from perturbation trajectories with valence
PTEST reference
PTPOST Averaging force field parameters over master and slaves
HCEST Resetting system coordinates to ab initio ref
HCEST Estimating force constants from Hessian cost for tasks HC_FC_DIAG HC_FC_
HCEST CROSS_ASS HC_FC_CROSS_ASA
HCEST Resetting system coordinates to ab initio ref
HCEST Estimating force constants from Hessian cost for tasks HC_FC_CROSS_ASS
HCEST HC_FC_CROSS_ASA HC_FC_CROSS_DSS HC_FC_CROSS_DSD HC_FC_CROSS_DAA HC_FC_
HCEST CROSS_DAD
BondHarm/C.C (kjmol/A^2 A)
fc = 3964 +- .000 rv = 1.394 +- .000
BondHarm/C.H (kjmol/A^2 A)
fc = 3330 +- .000 rv = 1.084 +- .000
BendAHarm/C.C.C (kjmol/rad^2 deg)
fc = 603.2 +- .000 rv = 119.9 +- .000
BendAHarm/C.C.H (kjmol/rad^2 deg)
fc = 320.9 +- .000 rv = 120.0 +- .000
TorsCheby2/C.C.C.C (kjmol au)
fc = 35.9 +- 0.0 sgn = -1
TorsCheby2/C.C.C.H (kjmol au)
fc = 31.2 +- 0.0 sgn = -1
TorsCheby2/H.C.C.H (kjmol au)
fc = 16.1 +- 0.0 sgn = -1
Oopdist/C.C.H.C (kjmol/A^2 A)
fc = 177.7 +- .000 rv = .0000 +- .000
Cross/C.C.C/b0a (kjmol/A A deg)
fc = 50.6 +- .0 rv0 = 1.39 +- 0.0 rv1 = 120 +- 0.0
Cross/C.C.C/b1a (kjmol/A A deg)
fc = 50.7 +- .0 rv0 = 1.39 +- 0.0 rv1 = 120 +- 0.0
Cross/C.C.C/bb (kjmol/A^2 A A)
fc = 525 +- .0 rv0 = 1.39 +- 0.0 rv1 = 1.39 +- 0.0
Cross/C.C.H/b0a (kjmol/A A deg)
fc = 114 +- .0 rv0 = 1.39 +- 0.0 rv1 = 119 +- 0.0
Cross/C.C.H/b1a (kjmol/A A deg)
fc = 138 +- .0 rv0 = 1.08 +- .00 rv1 = 119 +- 0.0
Cross/C.C.H/bb (kjmol/A^2 A A)
fc = 37.8 +- .0 rv0 = 1.39 +- 0.0 rv1 = 1.08 +- .00
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TIMING Initializing 0:00:00.101857
TIMING Equil Set RV 0:00:00.001024
TIMING PT Generate 0:00:05.651640
TIMING PT Estimate 0:00:00.059029
TIMING PT Post process 0:00:00.002952
TIMING HC Estimate FC 0:00:00.290666
__/\\\__________________________________________________________________/\\\____
\ \\\ \ \\\
\ \\\ End of file. Thanks for using QuickFF! Come back soon!! \ \\\
____\///__________________________________________________________________\///__
The logger will dump to following information to the screen (or a file if the
--logfile option was used):
- Machine environment:
- This section contains some general information about the machine environment on which you ran the quickff job. It contains the user name, machine info, starting time of job, Python version, Numpy version, Scipy versoin, Matplotlib version, current directory and the command used to envoke quickff.
- Routine sequence:
- Every task that is performed will be shown in the log output. For example,
the line
PTGEN Constructing trajectoriesindicates that QuickFF is busy constructing the perturbation trajectories. As a second example, the lineHCEST Estimating force constants from Hessian cost for tasks HC_FC_DIAG HC_FC_CROSSindicates QuickFF is estimating force constants by means of a least-square fit of the FF Hessian to the AI Hessian for each term for which the task HC_FC_DIAG and HC_FC_CROSS was selected.
- Final force field parameters:
- The final FF parameters are shown at the end of the log output.
- Timings:
- Finally, a summary of the timings for several steps during run of the program is provided.
The force field parameters were also written to the Yaff parameters file pars_yaff_noei.txt:
# BONDHARM
#---------
BONDHARM:UNIT K kjmol/A**2
BONDHARM:UNIT R0 A
BONDHARM:PARS C H 3.3306196160e+03 1.0844158109e+00
BONDHARM:PARS C C 3.9640168586e+03 1.3946450208e+00
# BENDAHARM
#----------
BENDAHARM:UNIT K kjmol/rad**2
BENDAHARM:UNIT THETA0 deg
BENDAHARM:PARS C C H 3.2099700494e+02 1.2000065910e+02
BENDAHARM:PARS C C C 6.0322228429e+02 1.1999902902e+02
# TORSION
#--------
TORSION:UNIT A kjmol
TORSION:UNIT PHI0 deg
TORSION:PARS C C C H 2 3.1236662655e+01 0.0000000000e+00
TORSION:PARS C C C C 2 3.5927820265e+01 0.0000000000e+00
TORSION:PARS H C C H 2 1.6155379141e+01 0.0000000000e+00
# OOPDIST
#--------
OOPDIST:UNIT K kjmol/A**2
OOPDIST:UNIT D0 A
OOPDIST:PARS C C H C 1.7779654600e+02 1.2582202304e-07
# Cross
#------
Cross:UNIT KSS kjmol/angstrom**2
Cross:UNIT KBS0 kjmol/(angstrom*rad)
Cross:UNIT KBS1 kjmol/(angstrom*rad)
Cross:UNIT R0 angstrom
Cross:UNIT R1 angstrom
Cross:UNIT THETA0 deg
Cross:PARS C C H 3.7836794362e+01 1.1489177691e+02 1.3847235970e+02 1.3946599160e+00 1.0844145084e+00 1.1999937065e+02
Cross:PARS C C C 5.2502130037e+02 5.0720412694e+01 5.0720412694e+01 1.3946599160e+00 1.3946599160e+00 1.2000556448e+02
Force field with electrostatics¶
- Generating Yaff input
First we need to generate the Yaff input file for the electrostatic contribution to the force field. This is done using the script qff_input-ei.py. For this tutorial, we will convert the charges given in the dataset
/chargesof the filegaussian_mbis.h5for the atoms in gaussian.fchk with atom types according to the level medium and use Gaussian charge distributions:qff-input-ei.py -v --ffatypes=low --gaussian gaussian.fchk gaussian_mbis.h5:charges pars_ei_mbisgauss.txt
This command dumped the following output to the screen, indicating wheter or not the atom types are well chosen from the point of view of electrostatics (see second remark in qff-input-ei.py):
Averaged Atomic Charges: ------------------------ H 0.131 +- 0.000 (N=6) C -0.132 +- 0.001 (N=6)Furthermore, the following Yaff parameter (pars_ei.txt) file was written:
#Fixed charges #--------------- FIXQ:UNIT Q0 e FIXQ:UNIT P e FIXQ:UNIT R angstrom FIXQ:SCALE 1 1.0 FIXQ:SCALE 2 1.0 FIXQ:SCALE 3 1.0 FIXQ:DIELECTRIC 1.0 # Atomic parameters # ---------------------------------------------------- # KEY label Q_0A R_A # ---------------------------------------------------- FIXQ:ATOM H 0.1314377213 0.7308000000 FIXQ:ATOM C -0.1315514661 1.1703000000
- Constructing the covalent contribution
Now, we generate a covalent force field on top of the previously derived electrostatic contribution using the qff.py script:
qff.py --ffatype=low --ei=pars_ei_mbisgauss.txt --suffix=_mbisgauss gaussian.fchk
The logging output for this job is:
________________/\\\_________/\\\\\\\\\\\\\\\__/\\\\\\\\\\\\\\\_________________ ______________/\\\\/\\\\_____\/\\\///////////__\/\\\///////////_________________ _____________/\\\//\////\\\___\/\\\_____________\/\\\___________________________ _____________/\\\______\//\\\__\/\\\\\\\\\\\_____\/\\\\\\\\\\\__________________ _____________\//\\\______/\\\___\/\\\///////______\/\\\///////__________________ _______________\///\\\\/\\\\/____\/\\\_____________\/\\\________________________ __________________\////\\\//______\/\\\_____________\/\\\_______________________ ______________________\///\\\\\\___\/\\\_____________\/\\\______________________ _________________________\//////____\///______________\///______________________ Welcom to QuickFF a Python package to quickly derive force fields from ab initio input data Written by Louis Vanduyfhuys(1)*, Steven Vandenbrande(1) and Toon Verstraelen(1) (1) Center for Molecular Modeling, Ghent University Belgium. * mailto: Louis.Vanduyfhuys@UGent.be USER louis MACHINE Linux T085 4.15.0-39-generic #42-Ubuntu SMP Tue Oct 23 MACHINE 15:48:01 UTC 2018 x86_64 TIME 2018-11-29 16:23:13.635020 QUICKFF VERSION 2.2.1 PYTHON VERSION 3.6.3 |Anaconda custom (64-bit)| (default, Oct 13 2017, PYTHON VERSION 12:02:49) [GCC 7.2.0] NUMPY VERSION 1.13.3 SCIPY VERSION 0.19.1 MATPLOTLIB VERSION 2.1.0 CURRENT DIR /home/louis/build/quickff/share/tutorials/benzene COMMAND LINE /home/louis/build/anaconda3/bin/qff.py --ffatypes=low COMMAND LINE --ei=pars_ei_mbisgauss.txt --suffix=_mbisgauss gaussian.fchk ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INIT Initializing system INIT Initializing Second order taylor reference for ai INIT Initializing Yaff force field reference for EI INIT Initializing program VAL Initializing valence force field VAL Added 12 bond terms VAL Added 18 bend terms (an)harmonic in the angle, 0 bend terms with VAL 1+cos(angle) potential and 0 bend terms with 1-cos(4*theta) potential. VAL Added 6 Harmonic and 0 SquareHarmonic out-of-plane distance terms EQSET Resetting system coordinates to ab initio ref EQSET Setting rest values to AI equilibrium values for tasks EQ_RV EQSET Averaging force field parameters over master and slaves PTGEN Resetting system coordinates to ab initio ref PTGEN Constructing trajectories PTEST Resetting system coordinates to ab initio ref PTEST Estimating FF parameters from perturbation trajectories PTPOST Averaging force field parameters over master and slaves VAL Resetting system coordinates to ab initio ref HCEST Resetting system coordinates to ab initio ref HCEST Estimating force constants from Hessian cost for tasks HC_FC_DIAG HC_FC_ HCEST CROSS_ASS HC_FC_CROSS_ASA PTEST Resetting system coordinates to ab initio ref PTEST Estimating FF parameters from perturbation trajectories with valence PTEST reference PTPOST Averaging force field parameters over master and slaves HCEST Resetting system coordinates to ab initio ref HCEST Estimating force constants from Hessian cost for tasks HC_FC_DIAG HC_FC_ HCEST CROSS_ASS HC_FC_CROSS_ASA HCEST Resetting system coordinates to ab initio ref HCEST Estimating force constants from Hessian cost for tasks HC_FC_CROSS_ASS HCEST HC_FC_CROSS_ASA HC_FC_CROSS_DSS HC_FC_CROSS_DSD HC_FC_CROSS_DAA HC_FC_ HCEST CROSS_DAD BondHarm/C.C (kjmol/A^2 A) fc = 3965 +- .000 rv = 1.394 +- .000 BondHarm/C.H (kjmol/A^2 A) fc = 3333 +- .000 rv = 1.087 +- .000 BendAHarm/C.C.C (kjmol/rad^2 deg) fc = 601.9 +- .000 rv = 119.9 +- .000 BendAHarm/C.C.H (kjmol/rad^2 deg) fc = 320.9 +- .000 rv = 120.1 +- .000 TorsCheby2/C.C.C.C (kjmol au) fc = 35.9 +- 0.0 sgn = -1 TorsCheby2/C.C.C.H (kjmol au) fc = 31.4 +- 0.0 sgn = -1 TorsCheby2/H.C.C.H (kjmol au) fc = 16.8 +- .00 sgn = -1 Oopdist/C.C.H.C (kjmol/A^2 A) fc = 166.2 +- .000 rv = .0000 +- .000 Cross/C.C.C/b0a (kjmol/A A deg) fc = 49.2 +- .0 rv0 = 1.39 +- 0.0 rv1 = 120 +- 0.0 Cross/C.C.C/b1a (kjmol/A A deg) fc = 49.3 +- .0 rv0 = 1.39 +- 0.0 rv1 = 120 +- 0.0 Cross/C.C.C/bb (kjmol/A^2 A A) fc = 525 +- .0 rv0 = 1.39 +- 0.0 rv1 = 1.39 +- 0.0 Cross/C.C.H/b0a (kjmol/A A deg) fc = 113 +- .0 rv0 = 1.39 +- .00 rv1 = 120 +- 0.0 Cross/C.C.H/b1a (kjmol/A A deg) fc = 139 +- .0 rv0 = 1.08 +- 0.0 rv1 = 120 +- 0.0 Cross/C.C.H/bb (kjmol/A^2 A A) fc = 40.3 +- .0 rv0 = 1.39 +- .00 rv1 = 1.08 +- 0.0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ TIMING Initializing 0:00:00.101447 TIMING Equil Set RV 0:00:00.001010 TIMING PT Generate 0:00:05.652945 TIMING PT Estimate 0:00:00.156266 TIMING PT Post process 0:00:00.003063 TIMING HC Estimate FC 0:00:00.311649 __/\\\__________________________________________________________________/\\\____ \ \\\ \ \\\ \ \\\ End of file. Thanks for using QuickFF! Come back soon!! \ \\\ ____\///__________________________________________________________________\///__An extra line appeared in the beginning of the log output, i.e.
INIT Initializing Yaff force field reference for EI. This indicates that an extra reference instance was created to represent the EI contribution to the force field. Furthermore, the covalent parameters are almost identical compared to the FF without electrostatics. This is indeed what we expect due to the charges being so small.The force field parameters were also written to the Yaff parameters file pars_yaff_mbisgauss.txt:
# BONDHARM #--------- BONDHARM:UNIT K kjmol/A**2 BONDHARM:UNIT R0 A BONDHARM:PARS C H 3.3330947115e+03 1.0876030857e+00 BONDHARM:PARS C C 3.9651506021e+03 1.3940024796e+00 # BENDAHARM #---------- BENDAHARM:UNIT K kjmol/rad**2 BENDAHARM:UNIT THETA0 deg BENDAHARM:PARS C C H 3.2097791046e+02 1.2015408283e+02 BENDAHARM:PARS C C C 6.0195035785e+02 1.1997128614e+02 # TORSION #-------- TORSION:UNIT A kjmol TORSION:UNIT PHI0 deg TORSION:PARS C C C H 2 3.1428590144e+01 0.0000000000e+00 TORSION:PARS C C C C 2 3.5982593481e+01 0.0000000000e+00 TORSION:PARS H C C H 2 1.6841711311e+01 0.0000000000e+00 # OOPDIST #-------- OOPDIST:UNIT K kjmol/A**2 OOPDIST:UNIT D0 A OOPDIST:PARS C C H C 1.6627879335e+02 1.4853054024e-07 # Cross #------ Cross:UNIT KSS kjmol/angstrom**2 Cross:UNIT KBS0 kjmol/(angstrom*rad) Cross:UNIT KBS1 kjmol/(angstrom*rad) Cross:UNIT R0 angstrom Cross:UNIT R1 angstrom Cross:UNIT THETA0 deg Cross:PARS C C H 4.0314253541e+01 1.1314940934e+02 1.3908936087e+02 1.3934733019e+00 1.0875970642e+00 1.2001585516e+02 Cross:PARS C C C 5.2510953701e+02 4.9272314421e+01 4.9272314421e+01 1.3934733019e+00 1.3934733019e+00 1.2003594739e+02
Tutorial 2 - Water¶
In the current tutorial, we will illustrate the usage of QuickFF as a library by constructing a basic force field for water. In Tutorial 3, we will cover the more advanced features of QuickFF.
For this tutorial we again assume to have the input files gaussian.fchk and
gaussian_mbis.h5. We will construct a script (which we will give the name
qff-derive.py) that will read the input files, convert them to data objects
QuickFF can read and implement a small program using QuickFF modules to derive
the force field. These files, together with the expected output, can be found
in share/tutorials/water.
First, we define the required settings using the Settings class:
#settings
with log.section('INIT', 2, timer='Initialization'):
settings = Settings(fn_yaff='pars_cov.txt', plot_traj='All', xyz_traj=True)
Second, we load the input files and convert them to a Yaff system instance and a QuickFF SecondOrderTaylor instance. Therefore, we use routines implemented in the fchk module of MolMod.
#load Gaussian Formatted Checkpoint file
fchk = FCHKFile('gaussian.fchk')
numbers = fchk.fields.get('Atomic numbers')
energy = fchk.fields.get('Total Energy')
coords = fchk.fields.get('Current cartesian coordinates').reshape([len(numbers), 3])
grad = fchk.fields.get('Cartesian Gradient').reshape([len(numbers), 3])
hess = fchk.get_hessian().reshape([len(numbers), 3, len(numbers), 3])
#Construct Yaff System file
system = System(numbers, coords)
system.detect_bonds()
system.set_standard_masses()
#Construct a QuickFF SecondOrderTaylor object containing the AI reference
with log.section('INIT', 2, timer='Initialization'):
ai = SecondOrderTaylor('ai', coords=coords, energy=energy, grad=grad, hess=hess)
At this instance, the system does not yet contain force field atom types. To define these, we could use the QuickFF feature guess_ffatypes. However, we can also use the ATSELECT language of Yaff:
#define atom types
rules = [
('H', '1 & =1%8'), #hydrogen atom with one oxygen neighbor
('O', '8 & =2%1'), #oxygen atom with two hydrogen neighbors
]
system.detect_ffatypes(rules)
The final input we need to parse, is the set of atomic charges. These charges
are read from the dataset charges in the HDF5 file guassian_mbis.h5.
#construct electrostatic force field from HE charges in gaussian_wpart.h5
with h5.File('gaussian_mbis.h5') as f:
charges = f['charges'][:]
scales = [1.0, 1.0, 1.0, 1.0]
with log.section('INIT', 2, timer='Initialization'):
ff_ei = get_ei_ff('EI', system, charges, scales, radii=None, average=True, pbc=[0,0,0])
scales is a list of 4 floats, representing the electrostatic scalings of 1-2,
1-3, 1-4 and 1-5 pairs. In this case, no scaling is applied to the electrostatic
interactions. radii=None indicates that point charges are used.
Average=True indicates that the set of charges are averaged over atoms of the same type.
Finally, pbc=[0,0,0] indicates that the system is not periodic. The returned
object ff_ei is an instance of the
YaffForceField class.
Now we have all required input, and we will write a new Program class, which will inherits from the BaseProgram class, and contain a sequence of instructions that will be executed consequently. In this case, first a diagonal force field will be constructed (i.e. without any cross terms), then we will add cross terms and fit the cross terms using a Hessian least squares fit while keeping the diagonal terms fixed:
class Program(BaseProgram):
def run(self):
with log.section('PROGRAM', 2):
#deriving diagonal force field
self.do_pt_generate()
self.do_pt_estimate()
self.average_pars()
self.do_hc_estimatefc(['HC_FC_DIAG'])
#adding and fitting cross terms
self.do_cross_init()
self.do_hc_estimatefc(['HC_FC_CROSS_ASS', 'HC_FC_CROSS_ASA'], logger_level=1)
#write output
self.make_output()
Finally, we can now setup the program and run it:
#initialize and run program
program = Program(system, ai, settings, ffrefs=[ff_ei] )
program.run()
To summarize, the entire script as well as the logger output is given below.
#!/usr/bin/env python
from yaff import System
from molmod.io.fchk import FCHKFile
from quickff.settings import Settings
from quickff.reference import SecondOrderTaylor, get_ei_ff
from quickff.program import BaseProgram
from quickff.log import log
import h5py as h5
#define class for deriving the force field
class Program(BaseProgram):
def run(self):
with log.section('PROGRAM', 2):
#deriving diagonal force field
self.do_pt_generate()
self.do_pt_estimate()
self.average_pars()
self.do_hc_estimatefc(['HC_FC_DIAG'])
#adding and fitting cross terms
self.do_cross_init()
self.do_hc_estimatefc(['HC_FC_CROSS_ASS', 'HC_FC_CROSS_ASA'], logger_level=1)
#write output
self.make_output()
#settings
with log.section('INIT', 2, timer='Initialization'):
settings = Settings(fn_yaff='pars_cov.txt', plot_traj='All', xyz_traj=True)
#load Gaussian Formatted Checkpoint file
fchk = FCHKFile('gaussian.fchk')
numbers = fchk.fields.get('Atomic numbers')
energy = fchk.fields.get('Total Energy')
coords = fchk.fields.get('Current cartesian coordinates').reshape([len(numbers), 3])
grad = fchk.fields.get('Cartesian Gradient').reshape([len(numbers), 3])
hess = fchk.get_hessian().reshape([len(numbers), 3, len(numbers), 3])
#Construct Yaff System file
system = System(numbers, coords)
system.detect_bonds()
system.set_standard_masses()
#Construct a QuickFF SecondOrderTaylor object containing the AI reference
with log.section('INIT', 2, timer='Initialization'):
ai = SecondOrderTaylor('ai', coords=coords, energy=energy, grad=grad, hess=hess)
#define atom types
rules = [
('H', '1 & =1%8'), #hydrogen atom with one oxygen neighbor
('O', '8 & =2%1'), #oxygen atom with two hydrogen neighbors
]
system.detect_ffatypes(rules)
#construct electrostatic force field from HE charges in gaussian_wpart.h5
with h5.File('gaussian_mbis.h5') as f:
charges = f['charges'][:]
scales = [1.0, 1.0, 1.0, 1.0]
with log.section('INIT', 2, timer='Initialization'):
ff_ei = get_ei_ff('EI', system, charges, scales, radii=None, average=True, pbc=[0,0,0])
#initialize and run program
program = Program(system, ai, settings, ffrefs=[ff_ei] )
program.run()
________________/\\\_________/\\\\\\\\\\\\\\\__/\\\\\\\\\\\\\\\_________________
______________/\\\\/\\\\_____\/\\\///////////__\/\\\///////////_________________
_____________/\\\//\////\\\___\/\\\_____________\/\\\___________________________
_____________/\\\______\//\\\__\/\\\\\\\\\\\_____\/\\\\\\\\\\\__________________
_____________\//\\\______/\\\___\/\\\///////______\/\\\///////__________________
_______________\///\\\\/\\\\/____\/\\\_____________\/\\\________________________
__________________\////\\\//______\/\\\_____________\/\\\_______________________
______________________\///\\\\\\___\/\\\_____________\/\\\______________________
_________________________\//////____\///______________\///______________________
Welcom to QuickFF
a Python package to quickly derive force fields from ab initio input data
Written by
Louis Vanduyfhuys(1)*, Steven Vandenbrande(1) and Toon Verstraelen(1)
(1) Center for Molecular Modeling, Ghent University Belgium.
* mailto: Louis.Vanduyfhuys@UGent.be
USER louis
MACHINE Linux T085 4.15.0-39-generic #42-Ubuntu SMP Tue Oct 23
MACHINE 15:48:01 UTC 2018 x86_64
TIME 2018-11-29 16:23:17.408563
QUICKFF VERSION 2.2.1
PYTHON VERSION 3.6.3 |Anaconda custom (64-bit)| (default, Oct 13 2017,
PYTHON VERSION 12:02:49) [GCC 7.2.0]
NUMPY VERSION 1.13.3
SCIPY VERSION 0.19.1
MATPLOTLIB VERSION 2.1.0
CURRENT DIR /home/louis/build/quickff/share/tutorials/water
COMMAND LINE qff-derive.py
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
INIT Initializing Second order taylor reference for ai
INIT Initializing Yaff force field reference for EI
INIT Initializing program
VAL Initializing valence force field
VAL Added 2 bond terms
VAL Added 1 bend terms (an)harmonic in the angle, 0 bend terms with
VAL 1+cos(angle) potential and 0 bend terms with 1-cos(4*theta) potential.
VAL Added 0 Harmonic and 0 SquareHarmonic out-of-plane distance terms
PTGEN Resetting system coordinates to ab initio ref
PTGEN Constructing trajectories
PTEST Resetting system coordinates to ab initio ref
PTEST Estimating FF parameters from perturbation trajectories
PROGRA Averaging force field parameters over master and slaves
HCEST Resetting system coordinates to ab initio ref
HCEST Estimating force constants from Hessian cost for tasks HC_FC_DIAG
VAL Resetting system coordinates to ab initio ref
HCEST Resetting system coordinates to ab initio ref
HCEST Estimating force constants from Hessian cost for tasks HC_FC_CROSS_ASS
HCEST HC_FC_CROSS_ASA
BondHarm/H.O (kjmol/A^2 A)
fc = 6247 +- .000 rv = 1.028 +- .000
BendAHarm/H.O.H (kjmol/rad^2 deg)
fc = 229.1 +- .000 rv = 91.53 +- .000
Cross/H.O.H/b0a (kjmol/A A deg)
fc = 78.5 +- .0 rv0 = 1.02 +- 0.0 rv1 = 91.5 +- 0.0
Cross/H.O.H/b1a (kjmol/A A deg)
fc = 78.5 +- .0 rv0 = 1.02 +- 0.0 rv1 = 91.5 +- 0.0
Cross/H.O.H/bb (kjmol/A^2 A A)
fc = -165 +- .0 rv0 = 1.02 +- 0.0 rv1 = 1.02 +- 0.0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TIMING Initializing 0:00:00.003614
TIMING Initialization 0:00:00.025545
TIMING PT Generate 0:00:00.181938
TIMING PT Estimate 0:00:00.005564
TIMING HC Estimate FC 0:00:00.016189
TIMING PT plot energy 0:00:00.503527
TIMING PT dump XYZ 0:00:00.000699
__/\\\__________________________________________________________________/\\\____
\ \\\ \ \\\
\ \\\ End of file. Thanks for using QuickFF! Come back soon!! \ \\\
____\///__________________________________________________________________\///__
Tutorial 3 - Biphenyl¶
In the current tutorial, we will illustrate advanced features of QuickFF by constructing a force field for biphenyl.
COMING SOON