5.4. Example: Pincer and Water Using topgen and Multiwfn
Tip
The sample input and output files can be found in testfiles/rigidmol/4-pincerwater.
In this section, we want to consider the cluster of a large molecular “pincer” and 13 water molecules. For water, we can just use tip4p.xyz as done before. The pincer is shown below:
98
pincer
O 0.97308000 4.02950000 0.02172960
O 5.22818000 4.29840000 0.55213100
O 5.45578000 2.88710000 -1.50102000
H 5.65508000 3.42840000 -0.65915200
O 1.63082000 1.39760000 -3.81621000
H 0.64178000 1.28790000 -3.63132000
O -0.85420000 1.13280000 -3.02372000
O -3.63480000 3.01800000 0.37844100
O 0.94203000 -4.12400000 -0.01164950
O 5.22756000 -4.12300000 -0.19236000
O 5.19975000 -2.75900000 1.90677000
H 5.49911000 -3.28700000 1.08741000
O 1.10372000 -1.43100000 3.83101000
H 0.13381000 -1.37300000 3.54797000
O -1.31340000 -1.31100000 2.81292000
O -3.65480000 -3.37200000 -0.81056300
N -4.55230000 1.09440000 -0.99695900
H -4.54060000 1.90820000 -0.33085600
N -4.76650000 -1.44200000 0.39144300
H -4.65210000 -2.26500000 -0.25439400
C 1.94344000 3.43050000 -0.80416400
C 3.32617000 3.54260000 -0.60306300
C 4.12217000 2.86750000 -1.60098000
C 3.56742000 2.16740000 -2.68707000
C 2.17005000 2.08920000 -2.79643000
C 1.35713000 2.72410000 -1.83511000
C -0.12610000 3.02720000 -1.80688000
C -1.12760000 1.90760000 -2.07830000
C -2.35500000 1.91710000 -1.32846000
C -2.54310000 2.90300000 -0.24390400
C -1.40410000 3.67780000 0.23983600
H -1.50880000 4.16770000 1.20107000
C -0.25120000 3.62020000 -0.44329800
C 3.98051000 4.29050000 0.46425400
C 3.18432000 5.08190000 1.47209000
H 3.87937000 5.54190000 2.17506000
H 2.48025000 4.43970000 2.00742000
H 2.58977000 5.85330000 0.97167700
C 4.46758000 1.52720000 -3.70897000
H 5.36536000 1.12230000 -3.23581000
H 4.80242000 2.26570000 -4.44793000
H 3.93826000 0.73680000 -4.24433000
C -0.45280000 4.13300000 -2.87362000
H -1.49930000 4.43910000 -2.78952000
H -0.26720000 3.72570000 -3.87041000
H 0.19468000 4.99770000 -2.70998000
C -3.38270000 0.95820000 -1.61680000
C -3.16450000 -0.17800000 -2.56830000
H -3.99600000 -0.88000000 -2.56379000
H -2.24950000 -0.71000000 -2.30253000
H -3.00750000 0.20520000 -3.58020000
C -5.75190000 0.27390000 -1.07838000
H -5.66800000 -0.41200000 -1.92466000
C -6.98160000 1.17200000 -1.28347000
H -7.04050000 1.88830000 -0.45262600
H -6.84510000 1.75430000 -2.20072000
C -8.26660000 0.33820000 -1.34828000
H -8.23510000 -0.30800000 -2.23605000
H -9.12970000 1.00130000 -1.46934000
C 1.79603000 -3.47200000 0.89376900
C 3.19591000 -3.51400000 0.81872300
C 3.86255000 -2.79400000 1.87848000
C 3.17612000 -2.10600000 2.89516000
C 1.76972000 -2.10300000 2.87705000
C 1.08107000 -2.79400000 1.86047000
C -0.37560000 -3.18400000 1.70563000
C -1.45940000 -2.12200000 1.86764000
C -2.61180000 -2.20900000 1.01188000
C -2.63740000 -3.19800000 -0.08655830
C -1.41550000 -3.90500000 -0.45472100
H -1.39180000 -4.38300000 -1.42651000
C -0.34400000 -3.78600000 0.33979400
C 3.97534000 -4.19500000 -0.20265100
C 3.31726000 -5.00200000 -1.29460000
H 4.09514000 -5.37900000 -1.95885000
H 2.60358000 -4.39700000 -1.86056000
H 2.75164000 -5.83600000 -0.86495000
C 3.94193000 -1.40900000 3.98735000
H 4.83366000 -0.92100000 3.58815000
H 4.27651000 -2.12900000 4.74422000
H 3.30842000 -0.67200000 4.48277000
C -0.73050000 -4.29800000 2.75283000
H -0.65690000 -3.87500000 3.75749000
H -0.02310000 -5.12500000 2.65752000
H -1.74630000 -4.66500000 2.58197000
C -3.69130000 -1.28100000 1.15965000
C -3.63830000 -0.14400000 2.13635000
H -4.52880000 0.47930000 2.09032000
H -2.76620000 0.47890000 1.93537000
H -3.50560000 -0.52900000 3.14983000
C -5.91880000 -0.57300000 0.21221300
H -6.00050000 0.11350000 1.05878000
C -7.19980000 -1.41500000 0.12954700
H -7.29880000 -1.99900000 1.05068000
H -7.10010000 -2.13200000 -0.69707000
C -8.42930000 -0.52700000 -0.09348310
H -9.32510000 -1.15200000 -0.17414200
H -8.57210000 0.12090000 0.78154800
5.4.1. Calculate Parameters
First, we will use topgen to generate most information:
$ topgen pincer.xyz
You will get the following output:
Analyze topology for the molecule from "pincer.xyz" ... done.
Exporting results to "pincer" +
.gjf: GJF file with bonding information.
-cycles.txt: Containing cycle texts that can be used for ABCluster:geom calculation.
-bonding.xyz: XYZ file with bonding information for ABCluster:geom calculation.
-rigid.xyz: XYZ file with bonding information for ABCluster:rigidmol calculation.
.psf: PSF file with topology informaition for ABCluster:geom/NAMD calculations.
In "pincer-rigid.xyz" and "pincer.psf":
"X" means unknown atomic type.
Total charge: 1.8360
In pincer.psf and pincer-rigid.xyz, the topology and Lennard-Jones parameters have been automatically generated. The atom typing is done using graph representation learning. Please refer to:
Tip
Zhang, J. Atom Typing Using Graph Representation Learning: How Do Models Learn Chemistry? J. Chem. Phys. 2022, 156, 204108.
However, the charge is 1.8360, which is wrong. To obtain reasonable charges, we recommend RESP charges, which can be calculated easily.
Step 1: Do a single point quantum chemistry calculation and save the wavefunction. For example, using Gaussian:
%nprocs=24
%mem=20GB
%chk=pincer.chk
#b3lyp/6-31g(d) geom=connectivity
title
0 1
[ coordinates ]
After calculation, run:
$ formchk pincer.chk pincer.fchk
pincer.fchk is an ASCII file and can be used in the following calculations.
Step 2: Get Multiwfn. Do as instructed below:
$ Multiwfn pincer.fchk
$ 7 // Select "Population analysis and calculation of atomic charges"
$ 18 // Select "Restrained ElectroStatic Potential (RESP) atomic charge"
$ 1 // Select "Start standard two-stage RESP fitting calculation"
Now it will calculate electrostatic potential at 48180 points for 1736 Gaussian basis functions. Thanks to the following algorithm, the calculation can be very fast:
Tip
Zhang, J.; Lu, T. Efficient Evaluation of Electrostatic Potential with Computerized Optimized Code. Phys. Chem. Chem. Phys. 2021, 23, 20323-20328.
After calculation, you can save the charges to pincer.chg. The last column is the RESP charges. Copy them to replace the charges in pincer.psf and pincer-rigid.xyz. The finished pincer-rigid.xyz is:
98
Generated by ABCluster
O 0.97308000 4.02950000 0.02172960
O 5.22818000 4.29840000 0.55213100
O 5.45578000 2.88710000 -1.50102000
H 5.65508000 3.42840000 -0.65915200
O 1.63082000 1.39760000 -3.81621000
H 0.64178000 1.28790000 -3.63132000
O -0.85420000 1.13280000 -3.02372000
O -3.63480000 3.01800000 0.37844100
O 0.94203000 -4.12400000 -0.01164950
O 5.22756000 -4.12300000 -0.19236000
O 5.19975000 -2.75900000 1.90677000
H 5.49911000 -3.28700000 1.08741000
O 1.10372000 -1.43100000 3.83101000
H 0.13381000 -1.37300000 3.54797000
O -1.31340000 -1.31100000 2.81292000
O -3.65480000 -3.37200000 -0.81056300
N -4.55230000 1.09440000 -0.99695900
H -4.54060000 1.90820000 -0.33085600
N -4.76650000 -1.44200000 0.39144300
H -4.65210000 -2.26500000 -0.25439400
C 1.94344000 3.43050000 -0.80416400
C 3.32617000 3.54260000 -0.60306300
C 4.12217000 2.86750000 -1.60098000
C 3.56742000 2.16740000 -2.68707000
C 2.17005000 2.08920000 -2.79643000
C 1.35713000 2.72410000 -1.83511000
C -0.12610000 3.02720000 -1.80688000
C -1.12760000 1.90760000 -2.07830000
C -2.35500000 1.91710000 -1.32846000
C -2.54310000 2.90300000 -0.24390400
C -1.40410000 3.67780000 0.23983600
H -1.50880000 4.16770000 1.20107000
C -0.25120000 3.62020000 -0.44329800
C 3.98051000 4.29050000 0.46425400
C 3.18432000 5.08190000 1.47209000
H 3.87937000 5.54190000 2.17506000
H 2.48025000 4.43970000 2.00742000
H 2.58977000 5.85330000 0.97167700
C 4.46758000 1.52720000 -3.70897000
H 5.36536000 1.12230000 -3.23581000
H 4.80242000 2.26570000 -4.44793000
H 3.93826000 0.73680000 -4.24433000
C -0.45280000 4.13300000 -2.87362000
H -1.49930000 4.43910000 -2.78952000
H -0.26720000 3.72570000 -3.87041000
H 0.19468000 4.99770000 -2.70998000
C -3.38270000 0.95820000 -1.61680000
C -3.16450000 -0.17800000 -2.56830000
H -3.99600000 -0.88000000 -2.56379000
H -2.24950000 -0.71000000 -2.30253000
H -3.00750000 0.20520000 -3.58020000
C -5.75190000 0.27390000 -1.07838000
H -5.66800000 -0.41200000 -1.92466000
C -6.98160000 1.17200000 -1.28347000
H -7.04050000 1.88830000 -0.45262600
H -6.84510000 1.75430000 -2.20072000
C -8.26660000 0.33820000 -1.34828000
H -8.23510000 -0.30800000 -2.23605000
H -9.12970000 1.00130000 -1.46934000
C 1.79603000 -3.47200000 0.89376900
C 3.19591000 -3.51400000 0.81872300
C 3.86255000 -2.79400000 1.87848000
C 3.17612000 -2.10600000 2.89516000
C 1.76972000 -2.10300000 2.87705000
C 1.08107000 -2.79400000 1.86047000
C -0.37560000 -3.18400000 1.70563000
C -1.45940000 -2.12200000 1.86764000
C -2.61180000 -2.20900000 1.01188000
C -2.63740000 -3.19800000 -0.08655830
C -1.41550000 -3.90500000 -0.45472100
H -1.39180000 -4.38300000 -1.42651000
C -0.34400000 -3.78600000 0.33979400
C 3.97534000 -4.19500000 -0.20265100
C 3.31726000 -5.00200000 -1.29460000
H 4.09514000 -5.37900000 -1.95885000
H 2.60358000 -4.39700000 -1.86056000
H 2.75164000 -5.83600000 -0.86495000
C 3.94193000 -1.40900000 3.98735000
H 4.83366000 -0.92100000 3.58815000
H 4.27651000 -2.12900000 4.74422000
H 3.30842000 -0.67200000 4.48277000
C -0.73050000 -4.29800000 2.75283000
H -0.65690000 -3.87500000 3.75749000
H -0.02310000 -5.12500000 2.65752000
H -1.74630000 -4.66500000 2.58197000
C -3.69130000 -1.28100000 1.15965000
C -3.63830000 -0.14400000 2.13635000
H -4.52880000 0.47930000 2.09032000
H -2.76620000 0.47890000 1.93537000
H -3.50560000 -0.52900000 3.14983000
C -5.91880000 -0.57300000 0.21221300
H -6.00050000 0.11350000 1.05878000
C -7.19980000 -1.41500000 0.12954700
H -7.29880000 -1.99900000 1.05068000
H -7.10010000 -2.13200000 -0.69707000
C -8.42930000 -0.52700000 -0.09348310
H -9.32510000 -1.15200000 -0.17414200
H -8.57210000 0.12090000 0.78154800
par_all36_prot q epsilon (kJ/mol) sigma (AA)
-0.1977931830 0.4184 2.9400 # OG3C51
-0.5733367878 0.5021 3.1181 # OG312
-0.5300303150 0.8037 3.1449 # OG311
0.4529362989 0.1925 0.4000 # HGP1
-0.4771591324 0.8037 3.1449 # OG311
0.3642538753 0.1925 0.4000 # HGP1
-0.4651447462 0.5021 3.0291 # OG2D1
-0.5666566754 0.5021 3.0291 # OG2D1
-0.1930103667 0.4184 2.9400 # OG3C51
-0.5726160092 0.5021 3.1181 # OG312
-0.5309161679 0.8037 3.1449 # OG311
0.4526673560 0.1925 0.4000 # HGP1
-0.4711202072 0.8037 3.1449 # OG311
0.3585367639 0.1925 0.4000 # HGP1
-0.4571705592 0.5021 3.0291 # OG2D1
-0.5623143971 0.5021 3.0291 # OG2D1
-0.5765288056 0.8368 3.2963 # NG2P1
0.3467211098 0.1925 0.4000 # HGP2
-0.5967315462 0.8368 3.2963 # NG2P1
0.3529095233 0.1925 0.4000 # HGP2
0.3879440331 0.2092 3.7418 # CG2R51
-0.5788065828 0.4142 3.3141 # CG2RC0
0.3842000880 0.1674 3.7418 # CG2R64
-0.0997233426 0.3766 3.3854 # CG2R62
0.3726833554 0.1674 3.7418 # CG2R64
-0.5441667574 0.2092 3.7418 # CG2R51
0.3735821875 0.1506 3.5814 # CG3C50
0.4739123225 0.4184 3.3854 # CG2R63
-0.7130252483 0.4142 3.3141 # CG2RC0
0.6954190395 0.4184 3.3854 # CG2R63
-0.5131154947 0.3766 3.3854 # CG2R62
0.2140815465 0.1925 1.9600 # HGR62
0.1989000377 0.2092 3.7418 # CG2R51
0.7609513260 0.4602 3.5636 # CG2O1
-0.4305219061 0.3264 3.6527 # CG331
0.1157827368 0.1004 2.3876 # HGA3
0.1157827368 0.1004 2.3876 # HGA3
0.1157827368 0.1004 2.3876 # HGA3
-0.1830077140 0.3264 3.6527 # CG331
0.0728685153 0.1004 2.3876 # HGA3
0.0728685153 0.1004 2.3876 # HGA3
0.0728685153 0.1004 2.3876 # HGA3
-0.3387068233 0.3264 3.6527 # CG331
0.0997306099 0.1004 2.3876 # HGA3
0.0997306099 0.1004 2.3876 # HGA3
0.0997306099 0.1004 2.3876 # HGA3
0.6233011950 0.4602 3.5636 # CG2N2
-0.4115073756 0.3264 3.6527 # CG331
0.1347767913 0.1004 2.3876 # HGA3
0.1347767913 0.1004 2.3876 # HGA3
0.1347767913 0.1004 2.3876 # HGA3
0.2434906893 0.3766 3.3854 # CG2R62
0.0278489446 0.1883 2.3876 # HGA1
-0.1453550327 0.3766 3.3854 # CG2R62
0.0609175806 0.1464 2.3876 # HGA2
0.0609175806 0.1464 2.3876 # HGA2
-0.0163527124 0.2803 3.5543 # CG2R71
0.0239680659 0.1464 2.3876 # HGA2
0.0239680659 0.1464 2.3876 # HGA2
0.3898962115 0.2092 3.7418 # CG2R51
-0.5779944634 0.4142 3.3141 # CG2RC0
0.3856230585 0.1674 3.7418 # CG2R64
-0.1005642641 0.3766 3.3854 # CG2R62
0.3698700817 0.1674 3.7418 # CG2R64
-0.5596427812 0.2092 3.7418 # CG2R51
0.4370629855 0.1506 3.5814 # CG3C50
0.4313162857 0.4184 3.3854 # CG2R63
-0.6910763065 0.4142 3.3141 # CG2RC0
0.6777322383 0.4184 3.3854 # CG2R63
-0.4926179545 0.3766 3.3854 # CG2R62
0.2114874971 0.1925 1.9600 # HGR62
0.1640672037 0.2092 3.7418 # CG2R51
0.7567282192 0.4602 3.5636 # CG2O1
-0.4286119953 0.3264 3.6527 # CG331
0.1151360088 0.1004 2.3876 # HGA3
0.1151360088 0.1004 2.3876 # HGA3
0.1151360088 0.1004 2.3876 # HGA3
-0.1807604294 0.3264 3.6527 # CG331
0.0722453737 0.1004 2.3876 # HGA3
0.0722453737 0.1004 2.3876 # HGA3
0.0722453737 0.1004 2.3876 # HGA3
-0.3417282498 0.3264 3.6527 # CG331
0.0988146354 0.1004 2.3876 # HGA3
0.0988146354 0.1004 2.3876 # HGA3
0.0988146354 0.1004 2.3876 # HGA3
0.6378338036 0.4602 3.5636 # CG2N2
-0.4424811015 0.3264 3.6527 # CG331
0.1442920859 0.1004 2.3876 # HGA3
0.1442920859 0.1004 2.3876 # HGA3
0.1442920859 0.1004 2.3876 # HGA3
0.2524345197 0.3766 3.3854 # CG2R62
0.0249875177 0.1883 2.3876 # HGA1
-0.1457390007 0.3766 3.3854 # CG2R62
0.0617790272 0.1464 2.3876 # HGA2
0.0617790272 0.1464 2.3876 # HGA2
-0.0227826499 0.2803 3.5543 # CG2R71
0.0255840761 0.1464 2.3876 # HGA2
0.0255840761 0.1464 2.3876 # HGA2
5.4.2. Global Optimization
Now, we can do the global optimization.
Step 1: Prepare an input file named pincer-wat.inp with the following content:
1pincer-wat.cluster # cluster file name
2100 # population size
3100 # maximal generations
43 # scout limit
510.0 # amplitude
6pincer-wat # save optimized configuration to .xyz and .gjf
730 # number of LMs to be saved
Step 2: Prepare the cluster file named pincer-wat.cluster with the following content:
12
2pincer-rigid.xyz 1
3tip4p.xyz 13
4* 20.0000
Step 3: Run the global optimization:
$ rigidmol pincer-wat.inp > pincer-wat.out
The global minimum is pincer-wat-OPT.xyz, which is shown below.