-
Notifications
You must be signed in to change notification settings - Fork 9
/
PyDuck.py
143 lines (122 loc) · 5.38 KB
/
PyDuck.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#!/usr/bin/env python
import os
import argparse
import pyscf
from pyscf import gto
import numpy as np
import subprocess
#Find the value of the environnement variable QUACK_ROOT. If not present we use the current repository
QuAcK_dir=os.environ.get('QUACK_ROOT','./')
#Create the argument parser object and gives a description of the script
parser = argparse.ArgumentParser(description='This script is the main script of QuAcK, it is used to run the calculation.\n If $QUACK_ROOT is not set, $QUACK_ROOT is replaces by the current directory.')
#Initialize all the options for the script
parser.add_argument('-b', '--basis', type=str, required=True, help='Name of the file containing the basis set in the $QUACK_ROOT/basis/ directory')
parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.')
parser.add_argument('-c', '--charge', type=int, default=0, help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0')
parser.add_argument('--cartesian', default=False, action='store_true', help='Add this option if you want to use cartesian basis functions.')
parser.add_argument('-fc', '--frozen_core', type=bool, default=False, help='Freeze core MOs. Default is false')
parser.add_argument('-m', '--multiplicity', type=int, default=1, help='Spin multiplicity. Default is 1 therefore singlet')
parser.add_argument('--working_dir', type=str, default=QuAcK_dir, help='Set a working directory to run the calculation.')
parser.add_argument('-x', '--xyz', type=str, required=True, help='Name of the file containing the nuclear coordinates in xyz format in the $QUACK_ROOT/mol/ directory without the .xyz extension')
#Parse the arguments
args = parser.parse_args()
input_basis=args.basis
unit=args.bohr
charge=args.charge
frozen_core=args.frozen_core
multiplicity=args.multiplicity
xyz=args.xyz + '.xyz'
cartesian=args.cartesian
working_dir=args.working_dir
#Read molecule
f = open(QuAcK_dir + '/mol/' + xyz,'r')
lines = f.read().splitlines()
nbAt = int(lines.pop(0))
lines.pop(0)
list_pos_atom = []
for line in lines:
tmp = line.split()
atom = tmp[0]
pos = (float(tmp[1]),float(tmp[2]),float(tmp[3]))
list_pos_atom.append([atom,pos])
f.close()
#Definition of the molecule
mol = gto.M(
atom = list_pos_atom,
basis = input_basis,
charge = charge,
spin = multiplicity - 1
)
#Fix the unit for the lengths
mol.unit=unit
#
mol.cart = cartesian
#Update mol object
mol.build()
#Accessing number of electrons
nelec=mol.nelec #Access the number of electrons
nalpha=nelec[0]
nbeta=nelec[1]
f = open(working_dir+'/input/molecule','w')
f.write('# nAt nEla nElb nCore nRyd\n')
f.write(str(mol.natm)+' '+str(nalpha)+' '+str(nbeta)+' '+str(0)+' '+str(0)+'\n')
f.write('# Znuc x y z\n')
for i in range(len(list_pos_atom)):
f.write(list_pos_atom[i][0]+' '+str(list_pos_atom[i][1][0])+' '+str(list_pos_atom[i][1][1])+' '+str(list_pos_atom[i][1][2])+'\n')
f.close()
#Compute nuclear energy and put it in a file
subprocess.call(['rm', working_dir + '/int/ENuc.dat'])
f = open(working_dir+'/int/ENuc.dat','w')
f.write(str(mol.energy_nuc()))
f.write(' ')
f.close()
#Compute 1e integrals
ovlp = mol.intor('int1e_ovlp')#Overlap matrix elements
v1e = mol.intor('int1e_nuc') #Nuclear repulsion matrix elements
t1e = mol.intor('int1e_kin') #Kinetic energy matrix elements
dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators
x,y,z = dipole[0],dipole[1],dipole[2]
norb = len(ovlp)
subprocess.call(['rm', working_dir + '/int/nBas.dat'])
f = open(working_dir+'/int/nBas.dat','w')
f.write(str(norb))
f.write(' ')
f.close()
def write_matrix_to_file(matrix,size,file,cutoff=1e-15):
f = open(file, 'a')
for i in range(size):
for j in range(i,size):
if abs(matrix[i][j]) > cutoff:
f.write(str(i+1)+' '+str(j+1)+' '+"{:.16E}".format(matrix[i][j]))
f.write('\n')
f.close()
#Write all 1 electron quantities in files
#Ov,Nuc,Kin,x,y,z
subprocess.call(['rm', working_dir + '/int/Ov.dat'])
write_matrix_to_file(ovlp,norb,working_dir+'/int/Ov.dat')
subprocess.call(['rm', working_dir + '/int/Nuc.dat'])
write_matrix_to_file(v1e,norb,working_dir+'/int/Nuc.dat')
subprocess.call(['rm', working_dir + '/int/Kin.dat'])
write_matrix_to_file(t1e,norb,working_dir+'/int/Kin.dat')
subprocess.call(['rm', working_dir + '/int/x.dat'])
write_matrix_to_file(x,norb,working_dir+'/int/x.dat')
subprocess.call(['rm', working_dir + '/int/y.dat'])
write_matrix_to_file(y,norb,working_dir+'/int/y.dat')
subprocess.call(['rm', working_dir + '/int/z.dat'])
write_matrix_to_file(z,norb,working_dir+'/int/z.dat')
#Write two-electron integrals
eri_ao = mol.intor('int2e')
def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
f = open(file, 'a')
for i in range(size):
for j in range(i,size):
for k in range(i,size):
for l in range(j,size):
if abs(tensor[i][k][j][l]) > cutoff:
f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
f.write('\n')
f.close()
subprocess.call(['rm', working_dir + '/int/ERI.dat'])
write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat')
#Execute the QuAcK fortran program
subprocess.call(QuAcK_dir+'/bin/QuAcK')