Growth of the Asymmetric Unit into a Crystal via Alchemy (GAUCHE)

The general procedure for calculating the organic crystal deposition/sublimation free energy using GAUCHE method is as follows:

1. Calculate the deposition free energy of the asymmetric unit.
2. Find the energy minimized crystal structure by optimizing the snapshots of the initial deposition simulation.
3. Calculate the free energy to add a harmonic restraint to each atom of the asymmetric unit between its current simulation coordinates and those of the energy minimized crystal.
4. Calculate the free energy change to transfer the restrained asymmetric unit molecules into vapor by alchemical annihilation.
5. Repeat the third and fourth steps using the unit cell instead of the asymmetric unit.

To do the deposition/sublimation free energy calculation, we need an input coordinate file in TINKER XYZ (or PDB) format. In this example, we are using acetanilide (ACANIL).

                    19  ACANIL.xyz
                    1    C   -1.62567502    4.09801640    2.43458664   401       2      11      12      13
                    2    C   -0.21564232    4.56925851    2.09474698   402       1       3       4
                    3    O    0.18378795    5.64372837    2.46911645   409       2
                    4    N    0.54676489    3.70436855    1.35104344   403       2       5      14
                    5    C    1.87236284    3.80268207    0.80739919   405       4       6      10
                    6    C    2.36589527    2.78623691   -0.02894495   404       5       7      15
                    7    C    3.66232337    2.84237728   -0.53485617   406       6       8      16
                    8    C    4.51553157    3.88888100   -0.21210739   407       7       9      17
                    9    C    4.02180870    4.91434631    0.58187883   406       8      10      18
                    10   C    2.71900496    4.89224328    1.07367960   404       5       9      19
                    11   H   -2.06287504    3.50784974    1.63545511   408       1
                    12   H   -2.29053235    4.93747119    2.60633250   408       1
                    13   H   -1.61366036    3.49096299    3.33493674   408       1
                    14   H    0.21515865    2.76325126    1.34137220   411       4
                    15   H    1.73971134    1.95334763   -0.27147787   410       6
                    16   H    4.00461342    2.04830814   -1.17749156   412       7
                    17   H    5.52996679    3.90951933   -0.56800284   413       8
                    18   H    4.65698146    5.74431910    0.83904803   412       9
                    19   H    2.38670936    5.70593100    1.68161589   410      10
                

The property file (like a TINKER keyword file) specifies the crystal information. Note that for most of the small organic molecules, AMOEBA parameters are not published thus need to be parameterized for each molecule (PolType Program developed by Pengyu Ren Research Group can be used for this purpose). Van der Waals cutoff is specified here, and the parameters and patch specifications point to acetanilide AMOEBA parameters ACANIL.patch.

                    spacegroup Pbca

                    a-axis    19.640
                    b-axis    9.483
                    c-axis    7.979
                    alpha     90.00
                    beta      90.00
                    gamma     90.00

                    vdw-cutoff  12.0

                    parameters ACANIL.patch
                    patch      ACANIL.patch
                

The command to calculate the deposition free energy of the asymmetric unit is as follows:

ffxc Thermodynamics -c 100 -a -l 0.0 -i stochastic -n 25000000 -t 298.15 -r 1.0 -s 1 -f 19 ACANIL.xyz

If experimental crystal structure is available, optimize the experimental structure to use for the remaining GAUCHE steps. Specify the RMS gradient convergence criteria using -e flag. Rename the output, i.e. energy minimized structure, as min.SG.xyz, or any naming convention that indicates that the coordinate file contains energy optimized asymmetric unit.

ffxc minimize -e 1e-4 ACANIL.xyz

To do the expansion of the asymmetric unit into a unit cell, some modifications to the properties files are necessary. For restrained crystal and restrained vapor, add "restrainterm true". For restrained vapor, also add "vdwterm false".

ffxc Thermodynamics -c 100 -a -l 2.0 -i stochastic -n 5000000 -t 298.15 -r 1.0 rminv.SG.xyz rmin.SG.xyz 

The command to calculate the free energy of adding a harmonic restraint to the asymmetric unit is as follows:

ffxc Thermodynamics -c 100 -a -l 2.0 -i stochastic -n 5000000 -t 298.15 -r 1.0 rmin.SG.xyz min.SG.xyz 

The command to calculate the free energy change to transfer the restrained asymmetric unit molecules into vapor by alchemical annihilation is as follows:

ffxc Thermodynamics -c 100 -a -l 2.0 -i stochastic -n 5000000 -t 298.15 -r 1.0 rminv.SG.xyz rmin.SG.xyz 

rminv.SG.xyz refers to the restrained vapor state of the asymmetric unit, and the rest of the commands are the same as the previous step.

Unit cell calculations are the same as the previous two steps. However, the coordinate file includes multiple copies of the molecule according to the number of symmetry operators for the spacegroup. The property file should also indicate the spacegroup of the molecule as P1. FFX command to generate a unit cell coordinate file is:

ffxc saveAsP1 min.SG.xyz