-
Notifications
You must be signed in to change notification settings - Fork 0
/
test.py
80 lines (59 loc) · 1.98 KB
/
test.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
import glob
from qmtools import BasisSet, Molecule, QMTools, Grid
import numpy
from ctypes import c_float
import time
# create a calculator object
calculator = QMTools()
basisset = BasisSet("cc-pvdz.bin")
# load the molecule,
# BCB: "./molecule_29766_0/"
# Benzene: "./molecule_241_0/"
# Hydrogen: "./molecule_783_0/"
# Water: "./molecule_962_0/"
folder = "./molecule_241_0/"
mol = Molecule(folder+"GEOM-B3LYP.xyz", folder+"D-CCSD.npy", basisset)
# the xyz must contain atomic coordinates in ANGSTROM
# - create a grid for the density
# step: grid step in ANGSTROM
# fat: 'fat' empty space around the molecule in ANGSTROM
step = 0.025
fat = 3.0
egrid = Grid.DensityGrid(mol, step, fat)
print(egrid)
density = calculator.ComputeDensity(mol, egrid)
calculator.WriteGrid_bin(mol, egrid, "density_"+str(step)+".bin")
#vgrid = Grid.MakeGrid([-8,-8,0], 0.1, [160, 160, 64])
#print()
#vgrid = Grid.DensityGrid(mol, 0.05, 3.0)
#print(vgrid)
#grid = Grid.TestGrid(mol, 0.1, 0)
#print(grid)
'''
grid.qube[1] = c_float(1)
grid.qube[grid.shape.x] = c_float(2)
grid.qube[grid.shape.x*grid.shape.y] = c_float(3)
print(grid._qube[1,0,0])
print(grid._qube[0,1,0])
print(grid._qube[0,0,1])
'''
# - parameters for potential calculations
# target_size: padded grid size,
# D: convolution constant
# Nd: number of convolutions
# target_size: size of padded grid, if size is kept egrid.shape.{xyz} no padding will be added
D = 0.01
Nd = 20
target_size = (egrid.shape.x, egrid.shape.y, egrid.shape.z)
print()
vgrid = egrid.CopyGrid()
calculator.ComputePotential_padded(mol, egrid, vgrid, target_size=target_size, spread="gauss")
calculator.WriteGrid_bin(mol, egrid, "potential_"+str(step)+".bin")
#density = calculator.ComputeDensity_subgrid(mol, egrid)
#calculator.WriteDensity(mol, egrid, "density_0.1_sg4.bin")
# t0 = time.time()
# hartree = calculator.ComputeHartree(mol, egrid, vgrid)
# print(f'Hartree solution time: {time.time()-t0:.4f}')
#calculator.WriteDensity(mol, vgrid, "hartree_exp.bin")
del mol
del calculator