Skip to content

JOYCE Tutorial. Part (I / IV)

Author: Anna Piras, Pablo M. Martinez

Contact: pablo.martinezm@uam.es

Keywords: Tutorial, FF, JOYCE

Welcome to the introductory JOYCE tutorial. Here, we will discuss the fundamentals of JOYCE and learn how to get set to build our own force fields (FF) for a simple flexible molecule. Part II will cover the fitting of stiff internal coordinates, part III will incorporate the flexible internal coordinates to conclude in part IV with the validation of the force field.

Our sample molecule and atom assignation

We will be developing a FF for the butane molecule. It is one of the simplest flexible system, useful as a first example.

Alt text

We should first decide how many chemically distinct atoms conform the molecule. This is usually assessed through symmetry arguments.

Alt text

The assigned names are:

  • C1 for the terminal carbon atoms
  • C2 for the inner carbon atoms
  • H1 for the hydrogen connected to the terminal carbon atoms
  • H2 for the hydrogen connected to the inner carbon atoms

QM data

JOYCE fits the force field on quantum-mechanical information. In the case of butane, and as in previous examples the hessian will be needed for the stiff internal coordinates. Additionally, some of the explored conformations by this molecule at room temperature cannot be described in the local harmonic approximation. As a result, the PES along each of the flexible rotations must be computed. Once the computations are finished, a formatted file (.fcc or .fchk) containing the relevant information must be created. In the case of flexible rotations, one of the obtained .fcc files is:

INFO
State file generated from file: butane.delta_000.fchk (format: fchk)

GEOM      UNITS=ANGS
    14
Geometry from butane.delta_000.fchk in xyz format (with filter: all)
C        1.47280266  -0.64773032  -0.00000000
C        0.78048231   0.72057489  -0.00000003
C       -0.78048229   0.72057491  -0.00000003
C       -1.47280267  -0.64773032  -0.00000001
H        2.56790140  -0.52118988   0.00000102
H        1.21293135  -1.24569633   0.88864110
H        1.21293314  -1.24569585  -0.88864197
H        1.13795746   1.28775722   0.87678508
H        1.13795752   1.28775730  -0.87678505
H       -1.13795740   1.28775720   0.87678512
H       -1.13795753   1.28775730  -0.87678502
H       -2.56790142  -0.52118982   0.00000006
H       -1.21293240  -1.24569622  -0.88864148
H       -1.21293213  -1.24569594   0.88864160

ENER      UNITS=AU
 -1.58449607E+02

GRAD      UNITS=AU
 -8.66364561E-06 -6.14388069E-05 -7.16997365E-09  2.43481416E-04 -3.17608234E-05
  5.84541573E-09 -2.43477969E-04 -3.17348998E-05 -2.16304644E-09  8.65668271E-06
 -6.14555048E-05 -2.71777979E-08 -3.44743928E-06  4.41933275E-06 -4.92080811E-09
  9.83928924E-06  1.27243327E-05  7.61248902E-07  9.87315919E-06  1.27071903E-05
 -7.62532178E-07 -4.25492041E-05  3.16781993E-05 -3.57586004E-05 -4.25541369E-05
  3.16695005E-05  3.57662046E-05  4.25557104E-05  3.16718192E-05 -3.57579382E-05
  4.25504127E-05  3.16634780E-05  3.57739901E-05  3.44740319E-06  4.42108740E-06
 -6.13821113E-09 -9.87041848E-06  1.27127150E-05 -7.52186621E-07 -9.84126012E-06
  1.27223798E-05  7.71538370E-07

Note how the structure is identical to those corresponding to the ground state, but in this case the hessian matrix is missing. This is because as the full scan is explored, we deviate from the LHA.

A computation for each structure must be calculated while spanning the whole rotation. This scan does not necessarily have to account for the 360 degrees if symmetry arguments can be applied. In this regard, a scan of 180 degrees suffices to capture all the PES variations.

Finally, through the atomic indexing, we can stabilish the equivalence between various rotations, i.e. all C2-C2-C1-H1 dihedrals are chemically indistinguisable. Thus, the same PES can be used to map all 6 dihedrals in the final QMD-FF.

Alt text

Note how, for the right plot, we could have forseen that the rotation would have a 3-fold degenerate minimum owing to the 3-fold rotation symmetry of the terminal C1 atom. In this regard, we could have performed a 60-degree scan, possibly with a smaller step between conformations to include more details of the PES. Finer scans require, naturally, more computational time to retrieve and the number of reasonable points is system and purpose-dependent.

Continue with the tutorial here to start the fitting procedure.