Hello psi4 forum
I would like to create a potential energy scan of two molecules with increasing distance using SAPT.
I have edited the z matrix to adjust the distance but at very small distances <1 Angstrom when visualising in Avogadro there was a difference between the z matrix and the 3d structure.
Should I simply vary distance using Z matrix or use a more elegant method with numpy and linear algebra?
Any advice is much appreciated!
The Z matrix should work.
Without a specific example with specific coordinates, there’s not much more we can say.
Thanks for you reply,
Have started with the below cartesian coordinates for dimer V, then converted them into z matrix using open babel, further below. Decreased the B20 variable that is H-O from 1.82 to 0.82 and converted back to cartesian, but finally measuring the distance in Avogadro gives 0.982 A
dimerV_original = psi4.geometry(“”"
nocom
noreorient
0 1
O 5.141300 -0.991500 0.002600
H 5.712100 -0.822000 -0.763400
H 5.702300 -0.820100 0.775400
“”")
#calculation
psi4.set_options({‘scf_type’: ‘df’,
‘freeze_core’: ‘true’})
psi4.energy(‘sapt0/jun-cc-pvdz’, molecule=dimerV_original)
dispV = psi4.variable(‘SSAPT0 DISP ENERGY’)
elstV = psi4.variable(‘SSAPT0 ELST ENERGY’)
exchV = psi4.variable(‘SSAPT0 EXCH ENERGY’)
indV = psi4.variable(‘SSAPT0 IND ENERGY’)
totV =psi4.variable(‘SSAPT0 TOTAL ENERGY’)
C
C 1 B1
C 1 B2 2 A2
C 2 B3 1 A3 3 D3
C 3 B4 1 A4 2 D4
C 5 B5 3 A5 1 D5
N 4 B6 2 A6 1 D6
C 7 B7 4 A7 2 D7
C 8 B8 7 A8 4 D8
H 9 B9 8 A9 7 D9
H 9 B10 8 A10 7 D10
H 9 B11 8 A11 7 D11
H 5 B12 3 A12 1 D12
H 6 B13 5 A13 3 D13
H 1 B14 2 A14 3 D14
O 3 B15 1 A15 2 D15
O 8 B16 7 A16 4 D16
H 2 B17 1 A17 3 D17
H 16 B18 3 A18 1 D18
H 7 B19 4 A19 2 D19
O 19 B20 16 A20 3 D20
H 21 B21 19 A21 16 D21
H 21 B22 19 A22 16 D22
Variables:
B1 1.39284
B2 1.40151
A2 119.97404
B3 1.40124
A3 121.13127
D3 359.99506
B4 1.40065
A4 118.91670
D4 0.00485
B5 1.39188
A5 121.21342
D5 359.99501
B6 1.41727
A6 117.57705
D6 180.00294
B7 1.37405
A7 129.05020
D7 180.00055
B8 1.52509
A8 114.43500
D8 179.99101
B9 1.09499
A9 114.42535
D9 0.01815
B10 1.09450
A10 108.67085
D10 238.21776
B11 1.09461
A11 108.67080
D11 121.81778
B12 1.08580
A12 118.54586
D12 180.00228
B13 1.08094
A13 120.84444
D13 179.99694
B14 1.08688
A14 120.10831
D14 180.00144
B15 1.36112
A15 122.93072
D15 180.00283
B16 1.22460
A16 124.38736
D16 359.99005
B17 1.08911
A17 119.09772
D17 180.00318
B18 0.98180
A18 109.74013
D18 0.00900
B19 1.01013
A19 114.66736
D19 359.99515
B20 1.82471
A20 173.38664
D20 181.83837
B21 0.97021
A21 107.38328
D21 302.24948
B22 0.97022
A22 107.13206
D22 54.53430
dimerV_minus1=psi4.geometry(“”"
nocom
noreorient
0 1
O -2.84203 -0.00266 0.34065
H -3.43763 0.76322 0.34182
H -3.42845 -0.77558 0.33612
“”" )
#calculation
psi4.set_options({‘scf_type’: ‘df’,
‘freeze_core’: ‘true’})
psi4.energy(‘sapt0/jun-cc-pvdz’, molecule=dimerV_minus1)
dispV_minus1 = psi4.variable(‘SSAPT0 DISP ENERGY’)
elstV_minus1 = psi4.variable(‘SSAPT0 ELST ENERGY’)
exchV_minus1 = psi4.variable(‘SSAPT0 EXCH ENERGY’)
indV_minus1 = psi4.variable(‘SSAPT0 IND ENERGY’)
totV_minus1 =psi4.variable(‘SSAPT0 TOTAL ENERGY’)
Is it psi4 or openbabel that is giving issues with the z-matrix?
When I use your z-matrix as a psi4 input, I can vary the B20 coordinate without issue.
I think its actually Avogadro, which I appreciate moves away from this discussion board.
Measuring the non covalent distance on a 3d molecule in Avogadro gives a different distance than any of the listed distances on the z matrix.
I split the z matrix into two and used a numpy array for the distance at B20 as above
A19 = 114.66736
D19 = 359.99515
--
0 1
O 19 B20 16 A20 3 D20
H 21 B21 19 A21 16 D21
H 21 B22 19 A22 16 D22
B20 = """ + str(distances[i]) + """
A20 = 173.38664
D20 = 181.83837
B21 = 0.97021
A21 = 107.38328
D21 = 302.24948
B22 = 0.97022
A22 = 107.13206
D22 = 54.53430
Hi all, thanks for your help
I have several paracetamol - water systems with different hydrogen bond arrrangement, I have the xyz and I want to run SAPT with varying hydrogen bond length.
When the oxygen of the water is bonded, I can use python to set a variable bond length ok using the z matrix. However when its the hydrogen on the water its giving me an issue.
Below is my original code, and then where I have attempted to change the z matrix for H23-O17 hydrogen bond water-carbonyl. The error that psi4 gives me is error finding anchor atom (I cant understand what the anchor atom wants to be?)
import psi4
import numpy as np
import matplotlib.pyplot as plt
psi4.core.clean()
psi4.set_memory(‘10 GB’)
distances = np.arange(1, 2.5, 1)
eelst = np.zeros((len(distances)))
eexch = np.zeros((len(distances)))
eind = np.zeros((len(distances)))
edisp = np.zeros((len(distances)))
esapt = np.zeros((len(distances)))
for i in range(len(distances)):
dimera = psi4.geometry(“”"
C
C 1 B1
C 1 B2 2 A2
C 2 B3 1 A3 3 D3
C 3 B4 1 A4 2 D4
C 5 B5 3 A5 1 D5
N 4 B6 2 A6 1 D6
C 7 B7 4 A7 2 D7
C 8 B8 7 A8 4 D8
H 9 B9 8 A9 7 D9
H 9 B10 8 A10 7 D10
H 9 B11 8 A11 7 D11
H 5 B12 3 A12 1 D12
H 6 B13 5 A13 3 D13
H 1 B14 2 A14 3 D14
O 3 B15 1 A15 2 D15
O 8 B16 7 A16 4 D16
H 2 B17 1 A17 3 D17
H 16 B18 3 A18 1 D18
H 7 B19 4 A19 2 D19
B1 = 1.39307
B2 = 1.39768
A2 = 119.78699
B3 = 1.40027
A3 = 120.95173
D3 = 0.38809
B4 = 1.39843
A4 = 119.38579
D4 = 359.54954
B5 = 1.39195
A5 = 121.01984
D5 = 0.03229
B6 = 1.42093
A6 = 117.14439
D6 = 178.42414
B7 = 1.36779
A7 = 129.17507
D7 = 158.88037
B8 = 1.52169
A8 = 114.90323
D8 = 181.64036
B9 = 1.09459
A9 = 114.26759
D9 = 5.58560
B10 = 1.09505
A10 = 108.71935
D10 = 243.97262
B11 = 1.09385
A11 = 108.57214
D11 = 127.55471
B12 = 1.08536
A12 = 118.83811
D12 = 180.69181
B13 = 1.08225
A13 = 119.00426
D13 = 180.88286
B14 = 1.08826
A14 = 119.80699
D14 = 179.71864
B15 = 1.36835
A15 = 122.92248
D15 = 179.74570
B16 = 1.22855
A16 = 124.22890
D16 = 2.01755
B17 = 1.08826
A17 = 119.17279
D17 = 180.63051
B18 = 0.96974
A18 = 109.05272
D18 = 0.16373
B19 = 1.01055
A19 = 114.63239
D19 = 344.49601
–
O 14 B20 6 A20 5 D20
H 21 B21 14 A21 6 D21
H 21 B22 14 A22 6 D22
##edited z matrix H23-O17 ###
O 23 B22 6 A20 5 D20
H 21 B21 14 A21 6 D21
H 17 B23 14 A22 6 D22
B20 = 2.38573
A20 = 172.25185
D20 = 325.93828
B21 = 0.96913
A21 = 93.88434
D21 = 327.57979
B22 = 0.97535
A22 = 63.75962
D22 = 225.20235
B23 = “”" + str(distances[i]) + “”"
""")
psi4.set_options({'scf_type': 'df',
'freeze_core': 'true',
})
#
# # save the plotsapt calculation
psi4.energy('sapt0/jun-cc-pvdz', molecule=dimera)
eelst[i] = psi4.variable('SAPT ELST ENERGY')
eexch[i] = psi4.variable('SAPT EXCH ENERGY')
eind[i] = psi4.variable('SAPT IND ENERGY')
edisp[i] = psi4.variable('SAPT DISP ENERGY')
esapt[i] = psi4.variable('SAPT TOTAL ENERGY')
psi4.core.clean()
plt.plot(distances, eelst, label=‘Electrostatics’)
plt.plot(distances, eexch, label=‘Exchange’)
plt.plot(distances, eind, label=‘Induction’)
plt.plot(distances, edisp, label=‘Dispersion’)
plt.plot(distances, esapt, label=‘SAPT Total’)
plt.xlabel(‘Distance (Angstrom)’)
plt.ylabel(‘Energy (kcal/mol)’)
plt.legend()
plt.savefig(‘dimera_sapt.png’)
data = np.column_stack((distances, eelst, eexch, eind, edisp, esapt))
np.savetxt(‘dimera_sapt.csv’, data, delimiter=‘,’,
header=‘Distance,Electrostatics,Exchange,Induction,Dispersion,SAPT Total’, comments=‘’)
plt.show()
Thanks!
It’s not clear to me what the original input is, and what the edit you made is. I understand that it has something to do with H23-O17, and I see the comment you made in the input file.
My best guess is that you’re saying you’re trying to change the first water into the second water. That won’t work because the oxygen is trying to refer to an atom that doesn’t exist yet. It’s atom 21, so it can only use atoms 1-20.
Also, please enclose your input files in backticks, per our advice on asking questions.
This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.