Quantum Cluster Growth¶
Contents
What is QCG?¶
Quantum Cluster Growth is a utility/driver extension of CREST for the xtb and xtb-IFF program. It is designed to automatically add explicit solvent molecules to a solute, to generate an ensemble of the resulting cluster, and to compute solvation free energies in an explicit fashion. Thereby, almost any solute–solvent combination can be used, including open-shell and charged solutes.
The key procedure is the automated addition of solvent molecules around the solute, mainly done by a screening of docking positions, using the xtb-IFF, and xtb geometry optimizations. A second step includes the ensemble generation, which can be done either with single a MD/MTD simulation or the NCI-MTD workflow of CREST. The solvation free energy can additionally be computed by generating a reference ensemble.
Note
Currently, a version of xtb-IFF and xtb has to be installed and to be accessible by calling xtb and xtbiff. Statically linked binaries of xtb-IFF can be found here and for xtb they can be found at the latest release page. A future release of xtb will make the use of xtb-IFF optional.
Automated Grow Algorithm and Ensemble Generation¶
The automated grow algorithm starts with input coordinates provided by the user. Also, additional information like the charge and the number of unpaired electrons of the solute can be given, as well as the number of solvent molecules that should be added. As performing an unbiased docking of the solvent would lead to solvent clustering, two wall potentials are applied.
Schematic depiction of the inner and outer wall potentials.¶
The inner one causes the solute to stay in the center, while the outer one prevents the solvent molecules from clustering and causes the solvent to surround the solute.
Having the prerequisites set, first an xtb-IFF docking of a solute and solvent molecule is performed under application of the wall potentials. Second, a geometry optimization is done with xtb. These two steps repeat until the user-defined number of solvents is added.
QCG Algorithm for the ensemble genereation.¶
Note
Without setting a number of solvents that should be added, QCG will add solvent molecules until the interaction energy converges. As this is often difficult, it might occur that the program will not stop adding solvents. Therefore, it is highly recommended to always define a number of solvent molecules that should be added.
The ensemble generation will also be performed with the wall potentials. It can be done with just an MD or MTD simulation and optimizing the snapshot geometries. Anyway, the NCI-MTD run-type is recommended and set as default because it explores the conformational space the most. It can be further enhanced by increasing the MTD time during this run-type or by decreasing the sampling frequency of snapshots.
Solvation Free Energy¶
The solvation free energy can also be computed with QCG in a supermolecular approach. To do so, a reference ensemble has to be generated first. By default, this is done with the Cut-Freeze-Fill (CFF) algorithm. It removes the solute from the highest populated clusters of the solute–solvent ensemble and fills the remaining cavity with solvents. Afterward, frequency calculations are performed for the solute–solvent and the reference ensembles, as well as the solute. These give rise to thermocorrections. After scaling the translational and rotational entropy and taking the conformational entropy into account, the free energies result. Substracting the free energy of the reference ensemble and the solute from the solute–solvent ensemble, yields finally the solvation free energy.
Supermolecular computation of the solvation free energy in QCG.¶
Note
The scaling factor for the translational and rotational entropy are empirically determined for all solvents. As this depends on the kind of solvent, it might be necessary to adjust this parameter. The printout always contains solvation free energies for different scaling factors, while the final result is given for the chosen scaling factor (default 0.75).
QCG Flags¶
The QCG extension of CREST is usually invoked via the command line. To do so, two files containing solute and solvent coordinates have to be provided that can be in any format supported by CREST. To activate QCG, use a call similar to
> crest [SOLUTE] -qcg [SOLVENT] [OPTIONS]
The general and technical options of CREST also apply for the QCG run-type ( CREST command line arguments). Take care to always set the number of cores with --T <INT>.
Algorithms¶
Algorithm |
Flag |
Description |
|---|---|---|
Grow [default] |
|
Use only the grow algorithm without ensemble generation |
Ensemble |
|
Use the grow algorithm with a subsequent ensemble generation |
Solvation Free Energy |
|
Generates a reference ensemble and computes the solvation free energy |
Growth¶
Flag |
Description |
|---|---|
|
Number of solvents that should be added |
|
Does not perform a GFN2-xTB preoptimization of the input structures |
|
Saves the tmp folders |
|
Use GFN1-xTB for geometry optimization during the Grow algorithm |
|
Use GFN2-xTB for geometry optimization during the Grow algorithm |
|
Use GFN-FF for geometry optimization during the Grow algorithm |
|
Use the same random number for the xTB-IFF runs |
|
Set the charge for the solute and creates .CHRG file |
|
Set the number of unpaired electrons for the solute and creates .UHF file |
|
Set the scaling factor for the outer wall potential. The default is 1.0 except for water |
|
Fix the solute during the grow process (automatically done for water as solvent) |
|
No fixing of the solute during the grow process (fixing is only applied for water as solvent file) |
Ensemble Generation¶
The defaults of the NCI-MTD run-type in QCG are the same as in the stand-alone use.
Flag |
Description |
|---|---|
|
Perform an ensemble generation with the NCI-MTD run-type |
|
Perform an ensemble generation with a single MTD simulation |
|
Perform an ensemble generation with a single MD simulation |
|
Use GFN1-xTB for M(T)D simulation and geometry optimizations during the ensemble generation |
|
Use GFN2-xTB for M(T)D simulation and geometry optimizations during the ensemble generation. |
|
Use GFN-FF for M(T)D simulation and geometry optimizations during the ensemble generation |
|
Set the M(T)D length |
|
Set the dumping frequency of the M(T)D simulations |
|
Set the M(T)D time step in fs |
|
Set the dumping frequency in which a new reference structure is taken for the bias potential |
|
Turn off the additional MDs on the lowest conformers after the MTD simulations in the NCI-MTD run-type |
|
Set the temperature for the additional MDs in the NCI-MTD run-type |
|
Maximum number of MTD restarts of the NCI-MTD run-type |
Solvation Free Energy¶
Flag |
Description |
|---|---|
|
Generate the reference ensemble with the same method the solute–solvent ensemble was generated instead of the CFF algorithm |
|
Compute only the pure solvation energy without frequency calculations |
|
Set the number of solute–solvent clusters to take further into account. The highest populated ones are considered. If not set, a number is determined according to the population distribution (maximal 10). |
|
Use GFN1-xTB for frequency calculations |
|
Use GFN2-xTB for frequency calculations |
|
Use GFN-FF for frequency calculations |
|
Set a scaling factor for the translational and rotational entropy. Only for printout. |
Example Applications¶
Building up Solvent Shells and using them for a Subsequent MD Simulation¶
Explicitly modeled solvent molecules can lead to different properties of the solute compared to implicit solvent models. Thus, it might be necessary to check on this, for example, for geometries and IR spectra.
Most stable gas-phase structure of Bacillaene with GFN-FF.¶
Let’s assume we have bacillaene and want to grow a cluster of 100 water molecules around it. Afterward, we perform an MD simulation to investigate the geometry in solution. To do so, we provide input coordinates of the solute bacillaene.xyz and of a water molecule water.xyz. As we do not have much time, GFN-FF is also used during the growth algorithm. We call CREST and activate the QCG algorithm with the following command.
crest bacillaene.xyz --qcg water.xyz --nsolv 100 --gfnff --T 12 --alpb water --nofix > crest.out
90
C -5.3127996594 -2.4157946011 -0.5090291244
C -6.6369198591 -2.2744765141 -0.3505867691
C -4.5376337067 -1.9989947690 -1.6511538708
C -7.4911082799 -1.4906802100 -1.3353795445
C -7.3417027536 -2.8130153570 0.8534300295
C -6.9820968905 -0.0260037238 -1.4231824713
C -3.2303223550 -1.6498120479 -1.6254976982
C -2.4648530707 -1.4782784301 -0.4345685359
C -1.1490119425 -1.1466912716 -0.3267799699
C -0.2345729980 -0.9625082034 -1.4997955626
C -0.6049768456 -0.9722141905 0.9890332154
C 0.6701604636 -0.6330961686 1.2870390935
C 1.1341096011 -0.4761647542 2.6269132552
C 2.3830259753 -0.1235992724 3.0110310022
C 3.4734843635 0.1679005598 2.1268262297
C 4.6931387856 0.5393278145 2.5267510375
C 5.8304304131 0.8677548326 1.6003866159
N 5.4283570547 0.8855111783 0.2104033711
C 5.3166273681 -0.1533114661 -0.6067536733
C 6.0134079054 -1.4661798013 -0.2726105716
O 4.7074640990 -0.0705419639 -1.6861449989
C 7.5144666461 -1.3152943375 -0.5941160582
O 5.4313708211 -2.5035661443 -1.0272412069
C -6.7935147347 0.4641208344 -0.0003944175
O -7.7403733383 0.5228053447 0.7742620645
N -5.5106869286 0.7636403711 0.3110847277
C -4.9463300150 1.1535464336 1.5274406520
C -3.6739854047 1.6059179412 1.5754484283
C -5.7878061510 1.0225637475 2.7522984886
C -2.7844664843 1.8243840709 0.4767015254
C -1.5364316121 2.3145689941 0.6208013650
C -0.6655220004 2.5888732490 -0.4843830298
C 0.5487325406 3.1343609346 -0.3631391694
C 8.3298335983 -2.5431938358 -0.1704211284
C 9.6891088069 -2.5276652589 -0.8678074623
C 8.5197655075 -2.5906603796 1.3443209467
O -8.8586444577 -1.4991160896 -1.0023350180
C 1.4399210852 3.4740558350 -1.5251049648
C 2.6756050844 2.6104389356 -1.3716483915
C 1.8111558123 4.9585806785 -1.5160495647
O 2.6206023923 1.5056225590 -2.0814961045
O 3.6041031244 2.8850885925 -0.6340527724
H -4.7625204190 -2.9056795986 0.2836333196
H -5.0509303359 -2.0042601899 -2.6049219243
H -7.4349047919 -1.9513459821 -2.3302603369
H -6.6760991929 -3.4329009729 1.4477280383
H -8.2048607459 -3.4019112683 0.5445011777
H -7.7116749441 -1.9955343308 1.4750390520
H -6.0619084015 0.0254594956 -1.9985287387
H -7.7560134130 0.5683074568 -1.9132455671
H -2.7489218428 -1.4394091975 -2.5706796840
H -3.0129717265 -1.5960544885 0.4909743144
H -0.7561184183 -1.0725573851 -2.4442035526
H 0.5638309381 -1.7039098370 -1.4666980415
H 0.2258577807 0.0242190429 -1.4708519245
H -1.2956296326 -1.1246557687 1.8098949728
H 1.3745898015 -0.4696821456 0.4862918222
H 0.3981454779 -0.6600516684 3.4004310551
H 2.5974468939 -0.0405121044 4.0685973172
H 3.2726215465 0.0815655312 1.0682418239
H 4.9262739703 0.6324234323 3.5775689496
H 6.6468639108 0.1562857796 1.7478921113
H 6.2119507967 1.8653636049 1.8407595221
H 4.8859889259 1.7105149385 -0.0683642472
H 5.8575220450 -1.7269366655 0.7774134480
H 7.6007536957 -1.1797242272 -1.6752148153
H 7.9082834268 -0.4205120198 -0.1104506371
H 5.0955504465 -2.0978322698 -1.8442489323
H -4.8389468734 0.6790182836 -0.4409166158
H -3.2785159256 1.8587339956 2.5484209564
H -5.2246095910 1.3178627663 3.6320617997
H -6.1254267820 -0.0065165747 2.8692478768
H -6.6790660370 1.6429351108 2.6701403518
H -3.1292138183 1.6120944762 -0.5273979730
H -1.1579767117 2.5402684410 1.6084450304
H -1.0491040613 2.3573854617 -1.4703055793
H 0.9505481083 3.3742128511 0.6120025660
H 7.7785462919 -3.4349097892 -0.4874864977
H 9.5694109024 -2.5429149313 -1.9488155636
H 10.2730806739 -3.3975457026 -0.5768929878
H 10.2450015890 -1.6328073712 -0.5947376585
H 9.0991604017 -3.4685164813 1.6200641190
H 7.5644928769 -2.6459407744 1.8604244387
H 9.0545665098 -1.7068730142 1.6874964687
H -8.9584252351 -0.9401749581 -0.2142721535
H 0.9490038267 3.1991238510 -2.4619332618
H 2.4126870178 5.2036544364 -2.3883354398
H 2.3855713372 5.1906177538 -0.6223783473
H 0.9080079216 5.5625777336 -1.5263116425
H 3.4359643823 0.9211814126 -1.9665264904
3
O -0.1918040235 1.3862489483 0.0047370042
H 0.7660977787 1.3911615443 -0.0141642652
H -0.4927337474 1.6150799341 -0.8756928250
==============================================
| |
| C R E S T |
| |
| Conformer-Rotamer Ensemble Sampling Tool |
| based on the GFN methods |
| P.Pracht, S.Grimme |
| Universitaet Bonn, MCTC |
==============================================
Version 2.11, Mon 19. Apr 11:43:20 CEST 2021
Using the xTB program. Compatible with xTB version 6.4.0
Cite work conducted with this code as
P. Pracht, F. Bohle, S. Grimme, PCCP, 2020, 22, 7169-7192.
and S. Grimme, JCTC, 2019, 15, 2847-2862.
with help from:
C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
C. Plett, P.Pracht, S. Spicher
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Command line input:
> crest bacillaene.xyz -qcg water.xyz -gfnff -nsolv 100 -T 20 -grow -alpb water -nofix
Solute-file: bacillaene.xyz
Solvent-file: water.xyz
-gfnff : Use of GFN-FF requested.
-T 20 (CPUs/Threads selected)
Use of GFN-FF for ensemble search requested.
Use of GFN-FF for frequency computation requested.
-mdtime 10 (MD length in ps)
--alpb water : implicit solvation
========================================
| ---------------- |
| Q C G |
| ---------------- |
| Quantum Cluster Growth |
| University of Bonn, MCTC |
========================================
S. Grimme, S. Spicher, C.Plett, unpublished.
=========================================
| quantum cluster growth: INPUT |
=========================================
QCG: Only Cluster Generation
input parameters
solute : bacillaene.xyz
charge : 0
uhf : 0
solvent : water.xyz
# of solvents to add : 100
# of cluster generated : 4
# of CPUs used : 20
Solvation model : water
xtb opt level : normal
System temperature [K] : 298.0
RRHO scaling factor : 0.75
Solute geometry
molecular radius (Bohr**1): 11.20
molecular area (Bohr**2): 2554.19
molecular volume (Bohr**3): 5887.65
Solvent geometry
molecular radius (Bohr**1): 3.88
molecular area (Bohr**2): 194.90
molecular volume (Bohr**3): 244.27
radius of solute : 18.06
radius of solvent : 6.25
=========================================
| Preoptimization |
=========================================
-------------------------
xTB Geometry Optimization
-------------------------
Geometry successfully optimized.
Generating LMOs for solute
Total Energy of solute: -127.4297277 Eh
-------------------------
xTB Geometry Optimization
-------------------------
Geometry successfully optimized.
Generating LMOs for solvent
Total energy of solvent: -5.0703134 Eh
________________________________________________________________________
__________________ Solute Cluster Generation _____________________
________________________________________________________________________
=========================================
| quantum cluster growth: GROW |
=========================================
Solute:
unit ellipsoid axis a,b,c : 0.428 0.289 0.283
Solvent:
unit ellipsoid axis a,b,c : 0.384 0.323 0.292
solvent anisotropy : 1.130
solute anisotropy : 1.197
roff inner wall : 6.997
solute max dist : 39.969
solvent max dist : 7.298
inner unit axis : 0.528 0.241 0.231
inner ellipsoid/Bohr: 36.893 16.810 16.173
outer ellipsoid/Bohr: 28.654 19.342 18.972
Size E /Eh De/kcal Detot/kcal Density Efix R av/act. Surface Opt
1 -132.538524 -24.15 -24.15 1.115 -13.744 0.0 0.0 6017.1 normal
2 -137.624631 -9.91 -34.06 1.117 -14.046 9.9 8.9 6188.1 normal
3 -142.708487 -8.50 -42.56 1.120 -14.343 9.4 8.8 6350.0 normal
4 -147.796548 -11.14 -53.69 1.119 -14.634 9.2 14.3 6536.7 normal
[...]
100 -636.352408 -7.47 ******* 1.198 -32.223 14.2 14.9 22281.9 normal
Final gfn2 optimization
Growth finished after 100 solvents added
Results can be found in grow directory
Energy list on file <qcg_energy.dat>
Interaction energy on file <qcg_conv.dat>
Growing process on <qcg_grow.xyz>
Final geometry after grow in <cluster.coord> and <cluster.xyz>
Potentials and geometry written in <cluster_cavity.coord> and <twopot_cavity.coord>
-----------------
Wall Time Summary
-----------------
Grow wall time : 0h :34m :32s
--------------------
Overall wall time : 0h :34m :32s
CREST terminated normally.
QCG automatically detects water as a solvent. This will cause the outer wall potential to be scaled to smaller sizes. Additionally, the solute will be fixed during the growth (also only in case of water). Bacillaene has a variety of different conformations and we want to consider the response of the intramolecular geometry upon addition of water. Thus, the --nofix flag was provided.
Note
Fixing the solute during the growth will increase the cluster quality, especially for water. It is only the default for water as solvent and can be switched off with --nofix. For other solvents, the solute can be fixed also by using the keyword --fixsolute.
Notice that we also choose the ALPB water model to get energies including an additional implicit solvation model around the cluster.
The resulting cluster can be found in the grow directory as cluster.xyz file. Additionally, each step is written to qcg_grow.xyz.
As we want to perform an MD simulation on this structure without dissociating the cluster, we also need the wall potentials found in the wall_potential file. This can be renamed to .xcontrol and directly used as an input for xtb to perform the constrained MD simulation.
Ensemble Generation and Microsolvation¶
Creating ensembles of generated clusters is important for various reasons. For example, the conformational space is explored during the used MD and MTD simulations so that new energy minima are usually found. Additionally, many problems require the weighting of different populated structures and the inclusion of the conformational entropy. As an example, a microsolvation approach is considered, but also large ensembles with multiple solvent shells can be generated similarly. As typically only a few solvents are added for this, the conformational space is rather small and it is possible to find relatively complete ensembles within a reasonable computational time. Now we want to add three water molecules to benzoic acid. For this, we again provide solute as well as solvent coordinates and call for the ensemble generation.
Tip
For larger clusters, the conformational space will increase. Therefore, the MTD time should be increased or the sampling frequency should be decreased. Using only one MD or MTD simulation will usually yield a much more incomplete ensemble.
crest benzoic_acid.xyz --qcg water.xyz --nsolv 3 --T 12 --ensemble --mdtime 50 --alpb water --wscal 1.0 --nofix > crest.out
15
H -5.151895 0.608937 0.184841
C -4.075803 0.560948 0.103703
C -3.304923 1.648961 0.482499
H -3.781062 2.542533 0.858155
C -1.927760 1.593624 0.380574
H -1.316613 2.433539 0.671921
C -1.315885 0.440886 -0.104813
C 0.159025 0.350784 -0.229059
O 0.718993 -0.633914 -0.685096
O 0.806733 1.411370 0.189344
C -2.093917 -0.650077 -0.484577
H -1.601704 -1.534740 -0.859582
C -3.469918 -0.588324 -0.379395
H -4.072688 -1.434587 -0.673879
H 1.807623 1.318950 0.057503
$coord
-0.36245704029697 2.61963060973384 0.00895163975603 O
1.44771485215846 2.62891406998886 -0.02676657950467 H
-0.93113174846333 3.05205846171614 -1.65481945499110 H
$end
==============================================
| |
| C R E S T |
| |
| Conformer-Rotamer Ensemble Sampling Tool |
| based on the GFN methods |
| P.Pracht, S.Grimme |
| Universitaet Bonn, MCTC |
==============================================
Version 2.11.1, Mon 16. Aug 09:59:32 CEST 2021
Using the xTB program. Compatible with xTB version 6.4.0
Cite work conducted with this code as
P. Pracht, F. Bohle, S. Grimme, PCCP, 2020, 22, 7169-7192.
and S. Grimme, JCTC, 2019, 15, 2847-2862.
with help from:
C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
C. Plett, P.Pracht, S. Spicher
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Command line input:
> crest benzoic_acid.xyz --qcg water.coord --nsolv 3 -T 12 -ensemble -mdtime 50 --alpb water --wscal 1.0 --nofix
Solute-file: benzoic_acid.xyz
Solvent-file: water.coord
-T 12 (CPUs/Threads selected)
-mdtime 50 (MD length in ps)
--alpb water : implicit solvation
========================================
| ---------------- |
| Q C G |
| ---------------- |
| Quantum Cluster Growth |
| University of Bonn, MCTC |
========================================
S. Grimme, S. Spicher, unpublished.
=========================================
| quantum cluster growth: INPUT |
=========================================
QCG: Cluster + Ensemble Generation
Ensemble generated via CREST
input parameters
solute : benzoic_acid.xyz
charge : 0
uhf : 0
solvent : water.coord
# of solvents to add : 3
Cluster generated that are above 10 % populated
# of CPUs used : 12
Solvation model : water
xtb opt level : normal
System temperature [K] : 298.0
RRHO scaling factor : 0.75
Solute geometry
molecular radius (Bohr**1): 6.57
molecular area (Bohr**2): 635.98
molecular volume (Bohr**3): 1188.36
Solvent geometry
molecular radius (Bohr**1): 3.88
molecular area (Bohr**2): 194.90
molecular volume (Bohr**3): 244.27
radius of solute : 10.59
radius of solvent : 6.25
=========================================
| Preoptimization |
=========================================
-------------------------
xTB Geometry Optimization
-------------------------
Geometry successfully optimized.
Generating LMOs for solute
Total Energy of solute: -26.1730317 Eh
-------------------------
xTB Geometry Optimization
-------------------------
Geometry successfully optimized.
Generating LMOs for solvent
Total energy of solvent: -5.0705444 Eh
________________________________________________________________________
__________________ Solute Cluster Generation _____________________
________________________________________________________________________
=========================================
| quantum cluster growth: GROW |
=========================================
Solute:
unit ellipsoid axis a,b,c : 0.408 0.306 0.286
Solvent:
unit ellipsoid axis a,b,c : 0.386 0.322 0.292
solvent anisotropy : 1.133
solute anisotropy : 1.169
roff inner wall : 1.388
solute max dist : 17.497
solvent max dist : 7.283
inner unit axis : 0.487 0.274 0.239
inner ellipsoid/Bohr: 14.890 8.363 7.292
outer ellipsoid/Bohr: 14.686 11.006 10.277
Size E /Eh De/kcal Detot/kcal Density Efix R av/act. Surface Opt
1 -31.277550 -21.32 -21.32 1.155 -7.372 0.0 0.0 1359.9 normal
2 -36.366081 -11.29 -32.61 1.143 -7.936 9.0 7.9 1551.0 normal
3 -41.458471 -13.71 -46.31 1.148 -8.463 9.1 10.0 1720.2 normal
Growth finished after 3 solvents added
Results can be found in grow directory
Energy list on file <qcg_energy.dat>
Interaction energy on file <qcg_conv.dat>
Growing process on <qcg_grow.xyz>
Final geometry after grow in <cluster.coord> and <cluster.xyz>
Potentials and geometry written in <cluster_cavity.coord> and <twopot_cavity.coord>
=========================================
| quantum cluster growth: ENSEMBLE |
=========================================
Method for ensemble search:--gff
Starting ensemble cluster generation by CREST routine
------------------------------------------------
Generating MTD length from a flexibility measure
------------------------------------------------
System flexiblity is set to 1.0 for NCI mode
flexibility measure : 1.000
t(MTD) / ps set by command line : 50.0
t(MTD) / ps : 50.0
Σ(t(MTD)) / ps : 600.0 (12 MTDs)
-------------------------------------
Starting a trial MTD to test settings
-------------------------------------
Estimated runtime for one MTD (50.0 ps) on a single thread: 1 min 15 sec
Estimated runtime for a batch of 12 MTDs on 12 threads: 1 min 15 sec
list of Vbias parameters applied:
$metadyn 0.00125 1.000
$metadyn 0.00083 1.000
$metadyn 0.00056 1.000
$metadyn 0.00125 0.667
$metadyn 0.00083 0.667
$metadyn 0.00056 0.667
$metadyn 0.00125 0.444
$metadyn 0.00083 0.444
$metadyn 0.00056 0.444
$metadyn 0.00125 0.296
$metadyn 0.00083 0.296
$metadyn 0.00056 0.296
*******************************************************************************************
** N E W I T E R A T I O N C Y C L E **
*******************************************************************************************
========================================
MTD Iteration 1
========================================
========================================
| Meta-MD (MTD) Sampling |
========================================
Starting Meta-MD 1 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0300
Vbias exp α /bohr⁻²: 1.00
Starting Meta-MD 2 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0200
Vbias exp α /bohr⁻²: 1.00
Starting Meta-MD 3 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0133
Vbias exp α /bohr⁻²: 1.00
Starting Meta-MD 4 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0300
Vbias exp α /bohr⁻²: 0.67
Starting Meta-MD 5 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0200
Vbias exp α /bohr⁻²: 0.67
Starting Meta-MD 6 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0133
Vbias exp α /bohr⁻²: 0.67
Starting Meta-MD 7 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0300
Vbias exp α /bohr⁻²: 0.44
Starting Meta-MD 8 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0200
Vbias exp α /bohr⁻²: 0.44
Starting Meta-MD 9 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0133
Vbias exp α /bohr⁻²: 0.44
Starting Meta-MD 12 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0133
Vbias exp α /bohr⁻²: 0.30
Starting Meta-MD 11 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0200
Vbias exp α /bohr⁻²: 0.30
Starting Meta-MD 10 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0300
Vbias exp α /bohr⁻²: 0.30
*Meta-MTD 8 finished*
*Meta-MTD 3 finished*
*Meta-MTD 5 finished*
*Meta-MTD 6 finished*
*Meta-MTD 4 finished*
*Meta-MTD 2 finished*
*Meta-MTD 1 finished*
*Meta-MTD 10 finished*
*Meta-MTD 12 finished*
*Meta-MTD 9 finished*
*Meta-MTD 11 finished*
*Meta-MTD 7 finished*
-----------------------
Multilevel Optimization
-----------------------
-------------------------
1. crude pre-optimization
-------------------------
Optimizing all 3000 structures from file "crest_rotamers_0.xyz" ...
1 [...] 3000
done.
input file name : crest_rotamers_1.xyz
output file name : crest_rotamers_2.xyz
reference state Etot : -4.01804455000000
3000 structures remain within 6.00 kcal/mol window
========================================
MTD Iteration 2
========================================
========================================
| Meta-MD (MTD) Sampling |
========================================
Starting Meta-MD 1 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0300
Vbias exp α /bohr⁻²: 1.00
[...]
Starting Meta-MD 9 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0133
Vbias exp α /bohr⁻²: 0.44
*Meta-MTD 3 finished*
*Meta-MTD 9 finished*
*Meta-MTD 7 finished*
*Meta-MTD 5 finished*
*Meta-MTD 1 finished*
*Meta-MTD 2 finished*
*Meta-MTD 6 finished*
*Meta-MTD 10 finished*
*Meta-MTD 4 finished*
*Meta-MTD 8 finished*
-----------------------
Multilevel Optimization
-----------------------
-------------------------
1. crude pre-optimization
-------------------------
Optimizing all 2500 structures from file "crest_rotamers_0.xyz" ...
1 [...] 2500
done.
input file name : crest_rotamers_1.xyz
output file name : crest_rotamers_2.xyz
reference state Etot : -4.01784771000000
2500 structures remain within 6.00 kcal/mol window
========================================
MTD Iterations done
========================================
Collecting ensmbles.
running RMSDs...
done.
E lowest : -4.01804
142 structures remain within 3.00 kcal/mol window
-----------------------------------------------
Additional regular MDs on lowest 3 conformer(s)
-----------------------------------------------
Starting MD 1 with the settings:
MD time /ps : 25.0
MD Temperature /K : 400.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
[...]
Starting MD 6 with the settings:
MD time /ps : 25.0
MD Temperature /K : 500.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
*MD 4 finished*
*MD 1 finished*
*MD 2 finished*
*MD 6 finished*
*MD 3 finished*
*MD 5 finished*
Appending file crest_rotamers_1.xyz with new structures
--------------------------------------------
Ensemble optimization with normal thresholds
--------------------------------------------
Optimizing all 892 structures from file "crest_rotamers_1.xyz" ...
1 [...] 892
done.
input file name : crest_rotamers_2.xyz
output file name : crest_rotamers_3.xyz
reference state Etot : -4.01951159000000
892 structures remain within 3.00 kcal/mol window
...............................................
A new lower conformer was found!
Improved by 0.00147 Eh or 0.92058kcal/mol
...............................................
*******************************************************************************************
** N E W I T E R A T I O N C Y C L E **
*******************************************************************************************
========================================
MTD Iteration 1
========================================
========================================
| Meta-MD (MTD) Sampling |
========================================
Starting Meta-MD 1 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0300
Vbias exp α /bohr⁻²: 1.00
[...]
Starting Meta-MD 9 with the settings:
MD time /ps : 50.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
dumpstep(Vbias)/ps : 1.0
Vbias factor k /Eh : 0.0133
Vbias exp α /bohr⁻²: 0.44
*Meta-MTD 3 finished*
*Meta-MTD 7 finished*
*Meta-MTD 5 finished*
*Meta-MTD 4 finished*
*Meta-MTD 1 finished*
*Meta-MTD 8 finished*
*Meta-MTD 2 finished*
*Meta-MTD 6 finished*
*Meta-MTD 9 finished*
*Meta-MTD 10 finished*
-----------------------
Multilevel Optimization
-----------------------
-------------------------
1. crude pre-optimization
-------------------------
Optimizing all 2500 structures from file "crest_rotamers_0.xyz" ...
1 [...] 2500
done.
input file name : crest_rotamers_1.xyz
output file name : crest_rotamers_2.xyz
reference state Etot : -4.01859099000000
2500 structures remain within 6.00 kcal/mol window
========================================
MTD Iterations done
========================================
Collecting ensmbles.
running RMSDs...
done.
E lowest : -4.01859
77 structures remain within 3.00 kcal/mol window
-----------------------------------------------
Additional regular MDs on lowest 3 conformer(s)
-----------------------------------------------
Starting MD 1 with the settings:
MD time /ps : 25.0
MD Temperature /K : 400.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
[...]
Starting MD 6 with the settings:
MD time /ps : 25.0
MD Temperature /K : 500.0
dt /fs : 1.5
dumpstep(trj) /fs : 200
*MD 5 finished*
*MD 6 finished*
*MD 4 finished*
*MD 2 finished*
*MD 1 finished*
*MD 3 finished*
Appending file crest_rotamers_1.xyz with new structures
--------------------------------------------
Ensemble optimization with normal thresholds
--------------------------------------------
Optimizing all 827 structures from file "crest_rotamers_1.xyz" ...
1 [...] 827
done.
input file name : crest_rotamers_2.xyz
output file name : crest_rotamers_3.xyz
reference state Etot : -4.01950240000000
827 structures remain within 3.00 kcal/mol window
================================================
| Final Geometry Optimization |
================================================
--------------------------------------------
Ensemble optimization with normal thresholds
--------------------------------------------
Optimizing all 827 structures from file "crest_rotamers_3.xyz" ...
1 [...] 827
done.
input file name : crest_rotamers_4.xyz
output file name : crest_rotamers_5.xyz
reference state Etot : -4.01950637000000
827 structures remain within 3.00 kcal/mol window
GFN2-xTB optimization
--------------------------------------------
Ensemble optimization with normal thresholds
--------------------------------------------
Optimizing all 827 structures from file "crest_rotamers_5.xyz" ...
1 [..] 827
done.
-------------------------------------------
Ensemble optimization with tight thresholds
-------------------------------------------
Optimizing all 827 structures from file "crest_rotamers_6.xyz" ...
1 [...] 827
done.
Single point computation with GBSA model
827 jobs to do.
done.
Cluster E /Eh Density Efix R av/act. Surface Opt
1 -41.458562 1.139 0.000 9.6 8.3 935.1 tight
[...]
827 -41.446970 1.121 0.000 6.5 7.9 926.7 tight
Conformers taken: 10
------------------------------------------------------------------------
------------------------------------------------------------------------
Boltz. averaged energy of final cluster:
G /Eh : -41.46409513
T*S /kcal : -1.364
Ensemble generation finished.
Results can be found in ensemble directory
Lowest energy conformer on file <crest_best.xyz>
List of full ensemble on file <full_ensemble.xyz>
List of used ensemble on file <final_ensemble.xyz>
Thermodynamical data on file <thermo_data>
Population of full ensemble on file <full_population.dat>
Population on file <population.dat>
-----------------
Wall Time Summary
-----------------
test MD wall time : 0h : 0m : 0s
MTD wall time : 0h : 0m :40s
multilevel OPT wall time : 0h : 2m :23s
MD wall time : 0h : 5m :56s
--------------------
Overall wall time : 0h : 9m : 8s
CREST terminated normally.
To make sure that we have a reasonable ensemble and energy minima, the MTD time was set to 50 ps. The ALPB solvent model was used to have a better energy ranking of the ensemble structures. It is only applied during final single-point computations. As the solvent is water, we used the --nofix flag so that the solute is not fixed during the Growth. Also, the scaling factor for the outer wall potential was set to 1.0.
Note
If water is used as a solvent input coordinate, special settings are applied during the cluster growth. The solute will be fixed in space and the outer wall potential will be adjusted to a smaller size. This safeguards reasonable structures during the growth process if complete solvent shells are desired. This causes the solvent to be added as close as possible to the origin, which is of course not always good for microsolvation. Therefore, the wall potential is set to larger values.
The result will be an ensemble, written to full_ensemble.xyz, the energetically lowest structure to crest_best.xyz, and a population of the clusters to full_population.dat.
Selected structures of the resulting ensemble with their relative energies in kcal/mol.¶
Constraining the Solute¶
Sometimes GFN2-xTB or GFN-FF geometry optimizations might distort DFT optimized structures. To prevent this, it is possible to constrain the solute geometry by providing an .xcontrol file.
Note
Constraining the solute with $fix will only work for the Growth, but not for the ensemble generations as it is not possible in xtb to fix atoms during M(T)Ds. Alternatively, the solute can be fixed with --fixsolute during the grow algorithm. This will not apply for the preoptimization.
Additionally, the charge can be set by the .CHRG file and the number of unpaired electrons by the .UHF file. The format has to be the same as for xtb or CREST. Alternatively, they can be evoked with the flags --chrg <INT> or --uhf <INT> that will create the respective files.
Note
If one of these files is present in the folder while calling CREST, they will be read automatically. So make sure that only the files are present that are required.
As an example, transition metal complexes tend to deform. One of them is shown in the following.
DFT optimized structure.¶
Now we want to constrain the ligands and to generate an ensemble. To do so, we provide the following .xcontrol file that constrains all bonds between the iron (16), the carbon atoms of the ring (atoms 3,4,6,7,8), and the CO ligands (atoms 17-22).
$constrain
atoms: 3,4,6-8,16,17-22
$end
Having a folder with this file named .xcontrol, the solute coordinates solute.xyz, and the solvent coordinates solvent.xyz present, the ensemble is now generated for example with
crest solute.xyz --qcg solvent.xyz --nsolv 6 --T 12 --ensemble --gbsa h2o --mdtime 50 --mddump 200 > crest.out
23
N 1.3802608000 -0.0318528000 0.0463356000
N -0.4099459000 -2.4279732000 -0.4426793000
C -0.8233287000 -1.1730691000 -0.1562094000
C 0.0282237000 -0.0329935000 0.0761397000
N -3.1128965000 1.1237167000 1.4357707000
C -0.8148438000 1.0531548000 0.5128170000
C -2.1758731000 -0.7690355000 0.1411386000
C -2.1257237000 0.5229799000 0.7675646000
H -0.4572325000 2.0269118000 0.8223111000
H -4.0464379000 0.7439366000 1.3737637000
H -3.0303702000 2.1068597000 1.6506659000
H -1.1169391000 -3.0848452000 -0.7424315000
H 0.4822106000 -2.5429436000 -0.9026828000
H 1.8456312000 -0.7155480000 -0.5345244000
H 1.8360609000 0.8698550000 0.0667192000
H -3.0359762000 -1.4262449000 0.1182750000
Fe -1.5405186000 0.5270943000 -1.3978036000
C -1.8580173000 -0.6232521000 -2.7083661000
C -2.9303295000 1.6215990000 -1.6740517000
C -0.3124958000 1.4435812000 -2.2886848000
O 0.5082387000 2.0393586000 -2.8321404000
O -2.0651590000 -1.4011695000 -3.5304385000
O -3.8189784000 2.3247399000 -1.8686288000
6
C -5.1936370000 1.8144547000 -0.0000255000
C -3.7637653000 1.6301290000 -0.0000193000
H -5.5336302000 2.0849167000 0.9967301000
H -5.6871686000 0.8944148000 -0.3037251000
H -5.4665453000 2.6066235000 -0.6930206000
N -2.6155337000 1.4820613000 0.0000603000
Two ensembles will result (a complete one as full_ensemble.xyz and a one with the highest populated clusters final_ensemble.xyz) with a solute structure that has the same Fe-C structure as the solute.xyz file.
A selected structure from the full_ensemble.xyz.¶
Warning
If the solute is constrained completely, the preoptimization of the solute will fail. Therefore, the preoptimization should be switched off with --nopreopt.
Solvation Free Energies¶
The solvation free energy can be computed for any solute-solvent combination. Again, only input geometries are required.
Let’s consider 1-pentanol in benzene. We provide the input coordinates pentanol.xyz and benzene.xyz. The following call will yield the solvation free energy.
crest pentanol.xyz --qcg benzene.coord --nsolv 25 --T 12 --gsolv --nclus 4 --fscal 0.65 --gbsa benzene > crest.out
18
C 1.1956067224 0.3810439760 0.2749699821
C 2.6993267582 0.1988054806 0.4398058012
C 3.3048739621 -0.5930561241 -0.7170407712
C 4.8104223372 -0.7810968097 -0.5556226318
C 5.4141490570 -1.5692919510 -1.7184358828
O 6.7925871816 -1.8090235819 -1.5591217978
H 0.7847307317 0.9460737154 1.1083269558
H 0.6967141904 -0.5847813396 0.2356859339
H 0.9746565928 0.9171147501 -0.6455272552
H 2.9009463902 -0.3234896023 1.3776341886
H 3.1790981167 1.1786519338 0.4972352731
H 3.1045811920 -0.0710584163 -1.6560239285
H 2.8240160617 -1.5719984791 -0.7751608993
H 5.0195701758 -1.3132022816 0.3749761863
H 5.2956555541 0.1982438034 -0.4961001357
H 5.2241018357 -1.0388860571 -2.6640512701
H 4.9541986817 -2.5581683607 -1.7828096756
H 7.2421309110 -0.9635515651 -1.4497510164
$coord
-13.87039726485370 2.82210248171752 -0.00000000000001 c
-12.01481772929320 4.66684409238559 0.00000000000005 c
-13.20059388975830 0.29275302437682 -0.00000000000006 c
-15.84108693251735 3.35633820529181 0.00000000000002 h
-9.48943548457286 3.98223577676367 -0.00000000000002 c
-12.53750229368448 6.64062893972117 0.00000000000005 h
-10.67521101562521 -0.39185500509649 0.00000000000005 c
-14.64860105080699 -1.14679576913850 -0.00000000000008 h
-8.81963207801176 1.45288625624259 0.00000000000002 c
-8.04142898121708 5.42178558198453 -0.00000000000009 h
-6.84894223900865 0.91865045782779 -0.00000000000001 h
-10.15253022931515 -2.36564093901184 0.00000000000007 h
$end
==============================================
| |
| C R E S T |
| |
| Conformer-Rotamer Ensemble Sampling Tool |
| based on the GFN methods |
| P.Pracht, S.Grimme |
| Universitaet Bonn, MCTC |
==============================================
Version 2.11, Mon 19. Apr 11:43:20 CEST 2021
Using the xTB program. Compatible with xTB version 6.4.0
Cite work conducted with this code as
P. Pracht, F. Bohle, S. Grimme, PCCP, 2020, 22, 7169-7192.
and S. Grimme, JCTC, 2019, 15, 2847-2862.
with help from:
C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
C. Plett, P.Pracht, S. Spicher
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Command line input:
> crest pentanol.xyz --qcg benzene.coord -nsolv 25 -T 12 -gsolv -nclus 4 -fscal 0.65 -gbsa benzene
Solute-file: pentanol.xyz
Solvent-file: benzene.coord
-T 12 (CPUs/Threads selected)
Use of GFN-FF for ensemble search requested.
Use of GFN-FF for frequency computation requested.
--gbsa benzene : implicit solvation
-mdtime 10 (MD length in ps)
========================================
| ---------------- |
| Q C G |
| ---------------- |
| Quantum Cluster Growth |
| University of Bonn, MCTC |
========================================
S. Grimme, S. Spicher, C.Plett, unpublished.
=========================================
| quantum cluster growth: INPUT |
=========================================
QCG: Calculation of delta G_solv
Ensemble generated via CREST
input parameters
solute : pentanol.xyz
charge : 0
uhf : 0
solvent : benzene.coord
# of solvents to add : 25
# of cluster generated : 4
# of CPUs used : 12
Solvation model : benzene
xtb opt level : normal
System temperature [K] : 298.0
RRHO scaling factor : 0.65
Solute geometry
molecular radius (Bohr**1): 6.52
molecular area (Bohr**2): 651.12
molecular volume (Bohr**3): 1159.46
Solvent geometry
molecular radius (Bohr**1): 6.06
molecular area (Bohr**2): 516.08
molecular volume (Bohr**3): 931.74
radius of solute : 10.51
radius of solvent : 9.77
=========================================
| Preoptimization |
=========================================
-------------------------
xTB Geometry Optimization
-------------------------
Geometry successfully optimized.
Generating LMOs for solute
Total Energy of solute: -20.8872146 Eh
-------------------------
xTB Geometry Optimization
-------------------------
Geometry successfully optimized.
Generating LMOs for solvent
Total energy of solvent: -15.8796407 Eh
________________________________________________________________________
__________________ Solute Cluster Generation _____________________
________________________________________________________________________
=========================================
| quantum cluster growth: GROW |
=========================================
[...]
=========================================
| quantum cluster growth: ENSEMBLE |
=========================================
[...]
________________________________________________________________________
_________________ Solvent Cluster Generation _____________________
________________________________________________________________________
Method for CFF: GFN2-xTB
=========================================
| quantum cluster growth: CFF |
=========================================
CUT-FREEZE-FILL Algorithm to generate reference solvent cluster
now adding solvents to fill cluster...
Size Cluster E /Eh De/kcal Detot/kcal Opt
------------------------------------------------------------------------
adding solvent is repulsive for cluster: 1
previous cluster taken...
26 2 -413.211286 -2.24 -2.24 tight
26 3 -413.210426 -3.44 -3.44 tight
26 4 -413.216554 -7.37 -7.37 tight
------------------------------------------------------------------------
volume filled
Starting optimizations + SP of structures
4 jobs to do.
done.
Cluster E /Eh Density Efix R av/act. Surface Opt
1 -397.332086 1.067 0.000 0.0 0.0 5955.5 tight
2 -397.318544 1.060 0.000 0.0 0.0 6373.5 tight
3 -397.317717 1.057 0.000 0.0 0.0 6591.9 tight
4 -397.323610 1.058 0.000 0.0 0.0 6433.2 tight
------------------------------------------------------------------------
------------------------------------------------------------------------
Boltz. averaged energy of final cluster:
G /Eh : -397.33208567
T*S /kcal : -0.001
Solvent cluster generation finished.
Results can be found in solvent_cluster directory
Structures on file <crest_ensemble.xyz>
Energies on file <cluster_energy.dat>
Population on file <population.dat>
=========================================
| quantum cluster growth: ESOLV |
| |
| -10.21 kcal/mol |
=========================================
=========================================
| Frequency evaluation |
=========================================
Method for CFF: GFN-FF
SOLUTE MOLECULE
Starting reoptimizations + Frequency computation of ensemble
1 jobs to do.
done.
SOLUTE CLUSTER
Starting reoptimizations + Frequency computation of ensemble
4 jobs to do.
done.
SOLVENT CLUSTER
Starting reoptimizations + Frequency computation of ensemble
4 jobs to do.
done.
Solute Gas properties
# H(T) SVIB SROT STRA G(T)
[kcal/mol] [ cal/mol/K ] [kcal/mol]
--------------------------------------------------------
108.59 20.37 27.20 39.32 82.68
Solute cluster properties
# H(T) SVIB SROT STRA G(T)
[kcal/mol] [ cal/mol/K ] [kcal/mol]
--------------------------------------------------------
1 1794.49 762.14 44.60 48.69 1539.44
2 1794.49 763.02 44.55 48.69 1539.19
3 1794.21 766.70 44.64 48.69 1537.80
4 1794.06 761.65 44.61 48.69 1539.16
Solvent cluster properties
# H(T) SVIB SROT STRA G(T)
[kcal/mol] [ cal/mol/K ] [kcal/mol]
--------------------------------------------------------
1 1683.97 718.09 44.20 48.55 1442.21
2 1683.88 721.14 44.26 48.55 1441.20
3 1683.17 730.73 44.55 48.55 1437.54
4 1683.61 722.85 44.32 48.55 1440.40
________________________________________________________________________
_________________________ Evaluation ____________________________
________________________________________________________________________
-----------------------------------------------------
Gsolv and Hsolv ref. state: [1 M gas/solution]
G_solv (incl.RRHO) = 3.65 kcal/mol
H_solv (incl.RRHO) = -7.64 kcal/mol
-----------------------------------------------------
-----------------------------------------------------
Gsolv and Hsolv ref. state: [1 M gas/solution]
G_solv (incl.RRHO) = 1.75 kcal/mol
H_solv (incl.RRHO) = -9.54 kcal/mol
-----------------------------------------------------
-----------------------------------------------------
Solvation free energies with scaled translational
and rotational degrees of freedom: Gsolv (scaling)
>> -16.93 (0.05) <<
>> -15.95 (0.10) <<
>> -14.97 (0.15) <<
>> -13.98 (0.20) <<
>> -13.00 (0.25) <<
>> -12.01 (0.30) <<
>> -11.03 (0.35) <<
>> -10.05 (0.40) <<
>> -9.06 (0.45) <<
>> -8.08 (0.50) <<
>> -7.10 (0.55) <<
>> -6.11 (0.60) <<
>> -5.13 (0.65) <<
>> -4.15 (0.70) <<
>> -3.16 (0.75) <<
>> -2.18 (0.80) <<
>> -1.20 (0.85) <<
>> -0.21 (0.90) <<
>> 0.77 (0.95) <<
>> 1.75 (1.00) <<
-----------------------------------------------------
==================================================
| Gsolv with SCALED RRHO contributions: 0.75 |
| [1 bar gas/ 1 M solution] |
| |
| G_solv (incl.RRHO)+dV(T)= -5.13 kcal/mol |
==================================================
-----------------
Wall Time Summary
-----------------
test MD wall time : 0h : 0m :13s
MTD wall time : 0h :39m : 3s
multilevel OPT wall time : 0h :46m :35s
MD wall time : 5h :23m :24s
CFF wall time : 0h :16m : 7s
Frequencies wall time : 0h : 1m :48s
--------------------
Overall wall time : 7h : 7m :12s
CREST terminated normally.
The call will cause the grow algorithm to start with subsequent ensemble generation. In addition to the solute-solvent ensemble also a pure solvent ensemble will be created from the solute-solvent ensemble. These reference clusters can be found in the folder solvent_ensemble.
The number of clusters that are considered for the solvation free energy were set to 4 with --nclus 4. This reduces the computational costs as only four reference clusters are computed and only 4 frequency calculations are performed per ensemble. Therefore, only the 4 energetically best clusters are written to final_ensemble.xyz.
Note
The reference ensemble is created per default with the CFF algorithm. This can be switched off by providing --nocff. Then, the second ensemble is generated similar to the solute-solvent cluster but without the solvent. However, this procedure is usually computationally more demanding and yields worse results. Therefore, it is recommended to always use the CFF algorithm.
The solvation free energy for pyridine in benzene at the given scaling factor will be printed. Additionally, the results for some other scaling factors are also shown.
Tip
The scaling factor of the translational and rotational entropy is empirically determined to be 0.75. As different solvents will quench these parts differently, also the scaling factor has to be adjusted. For example, 0.65 yielded better solvation free energies in case of benzene as solvent.
Hint
As finding an almost complete ensemble for a cluster containing many molecules is only feasible with a large computational effort, QCG suffers from a statistical error. As this is often no problem for many uses of the ensemble, the solvation free energy is in a way sensitive to the completeness. It is recommended to repeat the solvation free energy computation at least 5 times to reduce this uncertainty.
Restarting Computations¶
QCG has a restart functionality. If, for example, a grow directory is present from a previous QCG calculation, it will be read automatically. If the solute and solvent geometry, as well as the number of solvent molecules to add, match the cluster found in the directory, QCG will skip the grow algorithm and start with the ensemble search.
Warning
Be very careful. If the number of solvents to add is larger than the solvent molecules of the cluster in the grow directory, it will be deleted, as well as all other results files. Also, if directories are present that are a step further than the keyword, they will also be deleted. For example, using --grow in a directory where an ensemble folder is placed, will cause it to be deleted.
As an example, we again use benzoic acid from the previous example and add 30 water molecules with QCG
crest benzoic_acid.xyz --qcg water.xyz --nsolv 30 --T 12 --alpb water > crest.out
The resulting cluster looks good and we decide to generate an ensemble out of it. Using just the normal CREST NCI-MTD would yield a non-physical structure as the inner wall potential is missing. The hydrophilic moiety of the solute would thus move to the outer wall. Hence, we need to re-activate QCG. To do so, we have just to go into the directory we called QCG the last time and execute
crest benzoic_acid.xyz --qcg water.xyz --nsolv 30 --T 12 --alpb water --ensemble > crest.out
Tip
Important is that the input structures and the number of solvent given after --nsolv are matching the data in the grow directory. All other settings can be changed. Also, a solvation free energy computation can be done by substituting --ensemble by --gsolv.
QCG will read the cluster in the grow directory and start directly with the ensemble generation. After this, we can restart again with
crest benzoic_acid.xyz --qcg water.xyz --nsolv 30 --T 12 --alpb water --gsolv > crest.out
This will yield the solvation free energies.