Protonation Site Calculation
A guide to protonation/deprotonation sampling.
Table of contents
Protonation site screening
The protonation site screening is one of CREST’s original workflows. In the following, it is demonstrated for the alanineglycine molecule from Example 1.
Assuming the input coordinates are given as struc.xyz, then the screening procedure for the default settings of GFN2-xTB in the gasphase can be initiated via
crest input.tomlcrest struc.xyz --protonate
# This is a CREST input file
input = "struc.xyz"
runtype = "protonate"
threads = 4
[calculation]
[[calculation.level]]
method = "gfn2"
20
C 2.081440 0.615100 -0.508430
C 2.742230 1.824030 -1.200820
N 4.117790 1.799870 -1.190410
C 4.943570 2.827040 -1.822060
C 6.440080 2.569360 -1.637600
O 7.351600 3.252270 -2.069090
N 0.610100 0.695090 -0.538780
O 2.095560 2.724940 -1.739670
O 6.705220 1.463410 -0.897460
H 0.303080 1.426060 0.103770
H 0.338420 1.050680 -1.460480
C 2.488753 -0.593400 -1.198448
H 2.416500 0.557400 0.532050
H 4.614100 1.081980 -0.670550
H 4.699850 3.794460 -1.373720
H 4.722890 2.844690 -2.894180
H 7.687400 1.448620 -0.860340
H 2.029201 -1.457008 -0.719999
H 2.170233 -0.542411 -2.238576
H 3.572730 -0.688405 -1.154998
21
-33.9532959000000
C -2.2366791737 -0.5404887354 0.4121540226
C -1.0429992337 0.4322523796 0.1914454521
N 0.1273541796 0.1106459377 0.7023931033
C 1.3278207288 0.8751800087 0.4515985750
C 2.3614961324 0.0883604762 -0.3544268772
O 3.4019382844 0.5376718255 -0.7182378556
N -3.3679699719 0.2080825698 -0.1883769318
O -1.3010836947 1.4505236702 -0.4464138447
O 1.9536631778 -1.1652077323 -0.5882065573
C -2.0132179689 -1.8610502140 -0.3099228628
H -4.0811789893 0.5155580127 0.4739802483
H -2.7888404050 1.0835923166 -0.5679217435
H -2.4016755803 -0.6998826513 1.4799666105
H 0.2524997735 -0.7732256391 1.1756367513
H 1.7942383798 1.1763825594 1.3944326754
H 1.0646486020 1.7789046764 -0.1045882565
H 2.6272233518 -1.6496791341 -1.0942029622
H -2.8786309063 -2.5099583009 -0.1976369980
H -1.8273047541 -1.6929201908 -1.3685540258
H -1.1545333109 -2.3849925613 0.0991147410
H -3.8091519880 -0.2910958199 -0.9639833735
21
-33.9495757800000
C -2.2115960779 -0.4232578796 0.3354146276
C -0.9849427338 0.4641131361 0.1133708322
N 0.2286072373 0.0089398461 0.1068513485
C 1.3971435222 0.8377542711 -0.1024011044
C 2.6929687125 0.0343821360 0.0059880186
O 3.7734614371 0.5173993474 -0.1038635457
N -3.3398788397 0.4952076123 0.2242818787
O -1.2324295873 1.7060056707 -0.0721443475
O 2.4554161839 -1.2627297260 0.2339593620
C -2.2514600121 -1.5695854730 -0.6774397340
H -3.9570590462 0.4838042218 1.0250762361
H -3.8909606978 0.3287231337 -0.6099201196
H -2.1446928208 -0.8342484347 1.3494993462
H 0.4282588123 -0.9763775539 0.2587012566
H 1.4294852941 1.6512468588 0.6306086674
H 1.3655535266 1.3068129482 -1.0922428687
H 3.2887641565 -1.7590932864 0.2985989389
H -3.1563773476 -2.1533083300 -0.5276945700
H -2.2542621603 -1.1822765638 -1.6938410216
H -1.4031016048 -2.2366627382 -0.5583466487
H -2.2655920006 1.7583695831 -0.0105437454
21
-33.9075158100000
C -2.2251338634 -0.4294428722 0.3167057400
C -0.9673970868 0.3951686840 0.0040808989
N 0.1637246544 0.0496888934 0.6771325603
C 1.3462458434 0.8281808448 0.3749023518
C 2.4195420743 0.0456529289 -0.2450549336
O 3.6793896209 0.3562062589 -0.1680159273
N -3.4380171991 0.2195875500 -0.0687968484
O -0.9167601715 1.2971380917 -0.8032552399
O 2.1118427387 -1.0152233335 -0.9009340344
C -2.0887015559 -1.7600847916 -0.4348043739
H -3.6357897789 1.0602413023 0.4564787570
H -3.4734267738 0.4186111579 -1.0606659913
H -2.2610286904 -0.6128585416 1.3966539808
H 0.1829094047 -0.6824664346 1.3681234201
H 1.7079861224 1.4211977385 1.2234150460
H 0.8950301649 1.5380276922 -0.3873519278
H 2.8714577997 -1.4734236430 -1.3105708653
H -2.9880426576 -2.3445332030 -0.2657102207
H -1.9927834847 -1.5794005368 -1.5025590225
H -1.2274775786 -2.3257246098 -0.0910275619
H 3.9095477082 1.1696708650 0.3111684809
The production run in this example yields 3 possible structures of Ala-Gly+H+ at GFN2-xTB level in the gasphase. Two of the structures have a relative energy within 3 kcal/mol. The third structure has a much larger relative energy of almost 30 kcal/mol. These structures are written to the protonated.xyz ensemble file. Note, that the atom order of hydrogen atoms does change in between the structures and as such the ensemble does not formally match the ensemble file format, but due to the same atom type CREST can handle the file nonetheless.
Using other ions besides H+
It is possible to add other ions besides H+ with the --swel <symbol> command (short for “switch element”). To do so, simply specify the element symbol and the charge as <symbol> argument, for example
crest struc.xyz --protonate --swel ca2+
to add Ca2+ instead of H+. There is basically no limitation to which element/charge combination can be added with this command. However, adding polyatomic ions is currently not possible.
Deprotonation site screening
The deprotonation site screening is a very simple process. One-by-one, protons are removed from a given input structure and the resulting input geometries are optimized. In the following, it is demonstrated for the alanineglycine molecule from above .
Assuming the input coordinates are given as struc.xyz, then the screening procedure for the default settings of GFN2-xTB in the gasphase can be initiated via
crest input.tomlcrest struc.xyz --deprotonate
# This is a CREST input file
input = "struc.xyz"
runtype = "deprotonate"
threads = 4
[calculation]
[[calculation.level]]
method = "gfn2"
20
C 2.081440 0.615100 -0.508430
C 2.742230 1.824030 -1.200820
N 4.117790 1.799870 -1.190410
C 4.943570 2.827040 -1.822060
C 6.440080 2.569360 -1.637600
O 7.351600 3.252270 -2.069090
N 0.610100 0.695090 -0.538780
O 2.095560 2.724940 -1.739670
O 6.705220 1.463410 -0.897460
H 0.303080 1.426060 0.103770
H 0.338420 1.050680 -1.460480
C 2.488753 -0.593400 -1.198448
H 2.416500 0.557400 0.532050
H 4.614100 1.081980 -0.670550
H 4.699850 3.794460 -1.373720
H 4.722890 2.844690 -2.894180
H 7.687400 1.448620 -0.860340
H 2.029201 -1.457008 -0.719999
H 2.170233 -0.542411 -2.238576
H 3.572730 -0.688405 -1.154998
19
-33.5970120600000
C -2.1708593998 -0.5337349548 0.3174329502
C -1.1300518302 0.5631764406 -0.0109828000
N 0.1206150810 0.2013618108 0.2224174363
C 1.3321836188 0.9000800714 -0.0990976369
C 2.4989240716 -0.1817105892 0.0788604457
O 3.6394124385 0.1997854687 -0.1663904019
N -3.5145828210 0.0257215331 0.4393815808
O -1.5000115017 1.6529887352 -0.4378797093
O 2.0885706265 -1.2952519973 0.4542088317
C -2.1559906054 -1.5775218668 -0.7994319788
H -3.5368219047 0.6413357822 1.2486122040
H -3.6518357560 0.6345253879 -0.3671378904
H -1.9082787931 -1.0116891660 1.2677593853
H 0.4014134301 -0.7570338018 0.5027596031
H 1.5107597776 1.7382256168 0.5818246109
H 1.3263801551 1.2797019383 -1.1244436706
H -2.9077873850 -2.3398362648 -0.6051867349
H -2.3815551635 -1.1004724333 -1.7514818083
H -1.1764827299 -2.0438403226 -0.8636745012
19
-33.5584748600000
C -2.1927584335 -0.4616043612 0.3440786365
C -0.9552581293 0.3929737503 -0.0314541700
N 0.0746707491 0.1497799145 0.7312149580
C 1.2634058933 0.9097598712 0.4112450897
C 2.4176662088 0.0223878739 0.0810903851
O 3.5192313290 0.0557773687 0.5950471942
N -3.4335471612 0.2433576147 0.0009454232
O -1.0793986160 1.2000279714 -0.9808987844
O 2.1563918940 -0.8800339117 -0.8845041904
C -2.1334963814 -1.7717852649 -0.4380978387
H -3.5464154298 1.0130709993 0.6585739972
H -3.2245641905 0.6940212025 -0.8953163391
H -2.1778559810 -0.6749250459 1.4174900496
H 1.5609611501 1.5455410066 1.2583048133
H 1.0643516601 1.5513773159 -0.4696834509
H 2.9336605171 -1.4482855323 -0.9998771010
H -3.0010004151 -2.3943363262 -0.2161852163
H -2.1207001621 -1.5550657049 -1.5047051136
H -1.2222575076 -2.3037500664 -0.1775085954
19
-33.5580574900000
C -2.1690752402 -0.5896911445 0.3007278255
C -1.0770071421 0.4696696924 0.0142832995
N 0.1623498042 0.0267699146 0.1855561818
C 1.3802255595 0.7017906623 -0.0182948372
C 2.5626404067 0.0611685908 0.2016639002
O 3.7297155720 0.4777116712 0.0801592636
N -3.4599663836 0.0497353208 0.5540210958
O -1.4116252199 1.6071148069 -0.3495343490
O 2.4182039037 -1.3021064678 0.6393570478
C -2.3049547011 -1.5169584430 -0.9091156655
H -3.4010727265 0.5478312871 1.4393397960
H -3.5485107853 0.7806781467 -0.1540026061
H -1.8977951026 -1.1816934881 1.1843481101
H 0.3058623525 -0.9263106819 0.4897142094
H 1.3253072133 1.7195073133 -0.3521321533
H 3.3270672171 -1.6049746565 0.7469409798
H -3.0981100285 -2.2434691439 -0.7396358968
H -2.5528476412 -0.9284333481 -1.7899106278
H -1.3689397936 -2.0398200288 -1.0912306054
The deprotonation yields 3 possible structures of Ala-Gly- at GFN2-xTB level in the gasphase. Only one of these structures is favorable, while the other two isomers have relative energies at about 24 kcal/mol. These structures are written to the deprotonated.xyz ensemble file. Note, that like for the protonation procedure, the atom order of hydrogen atoms does change in between the structures and as such the ensemble does not formally match the ensemble file format.
Tautomerism screening
The third application belonging to the protonation/deprotonation procedures described in this example is the screening for prototropic tautomers. This can be done simply by consecutive execution of protonation and deprotonation site screening for a molecule. In the following, it is demonstrated for the alanineglycine molecule from above .
Assuming the input coordinates are given as struc.xyz, then the screening procedure for the default settings of GFN2-xTB can be initiated as before. However, since we expect zwitter ions as the possible tautomers of Ala-Gly, we also use ALPB implicit solvation (for water) in this example.
crest input.tomlcrest struc.xyz --tautomerize --alpb water
# This is a CREST input file
input = "struc.xyz"
runtype = "tautomerize"
threads = 4
[calculation]
[[calculation.level]]
method = "gfn2"
alpb = "water"
20
C 2.081440 0.615100 -0.508430
C 2.742230 1.824030 -1.200820
N 4.117790 1.799870 -1.190410
C 4.943570 2.827040 -1.822060
C 6.440080 2.569360 -1.637600
O 7.351600 3.252270 -2.069090
N 0.610100 0.695090 -0.538780
O 2.095560 2.724940 -1.739670
O 6.705220 1.463410 -0.897460
H 0.303080 1.426060 0.103770
H 0.338420 1.050680 -1.460480
C 2.488753 -0.593400 -1.198448
H 2.416500 0.557400 0.532050
H 4.614100 1.081980 -0.670550
H 4.699850 3.794460 -1.373720
H 4.722890 2.844690 -2.894180
H 7.687400 1.448620 -0.860340
H 2.029201 -1.457008 -0.719999
H 2.170233 -0.542411 -2.238576
H 3.572730 -0.688405 -1.154998
In general, a much larger variety of structures is obtained than for the standalone protonation or deprotonation procedure. For Ala-Gly, as expected the zwitter ion is the most stable tautomer in implicit solvation. Other tautomers, such as the (S)-isomer to the input (R)-structure are also generated. However, since many chemical changes can occur in the protonation/deprotonation sequence, also many artifacts at the respective level of theory are created. When evaluating the output file (tautomers.xyz), the user must decide by his/her chemical intuition whether to consider a generated tautomer for further investigation or not.
Note, that the protonation/deprotonation sequence is performed twice for the tautomer screening. This is because a single sequence swaps one proton position relative to the input structure. Swapping more than one proton position squentially can therefore be achieved by multiple executions of the protonation/deprotonation.