-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.py
123 lines (96 loc) · 5.48 KB
/
main.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
import sys, os, subprocess
from sys import argv
from math import sqrt, floor, ceil
# import all user-defined functions for working with water files
import clusters
import clusterdata
import recursive
# =============================================================
# ------------- the main program ------------------------------
# =============================================================
# read command line arguments
script, infile = argv
clusterdata.Rcutoff=float(clusterdata.Rcutoff)
#filename = os.path.basename(infile)
clusterdata.snapshotdir = os.path.dirname(infile)
# read the box size from the cp2k input template
# the size of this box determines the shift in gas_phase calculations
clusterdata.abc_gasphase[:]=clusters.get_abc_from_template(clusterdata.cp2kfname+".inp")
# determine the largest cluster we have to do
isize=0
for icluster in clusterdata.docluster:
isize += 1
if (icluster==1):
clusterdata.largest_cluster = isize
if (clusterdata.largest_cluster < 1):
exit(3)
# check if mode_epsilon_file array is the same length as docluster
if (len(clusterdata.docluster) != len(clusterdata.mode_epsilon_file) ):
print "Lenght of docluster and mode_epsilon_file must be equal"
exit(3)
clusterdata.ncells = [0,0,0] # how many periodic cells in each direction we should consider
for idim in range(0,3):
clusterdata.ncells[idim] = int(ceil( ((clusterdata.largest_cluster-1)*clusterdata.Rcutoff)/clusterdata.abc[idim] ))
# it is important that Rcutoff is such that the largest clusters is less than half of a cell
# otherwise spurious convergence efffect will be observed
if ( (clusterdata.largest_cluster-1)*clusterdata.Rcutoff > clusterdata.abc[idim]/2.0 ):
print "Rcutoff is too large for this system. Rcutoff must be smaller than %10.5f" % (clusterdata.abc[idim]/(2.0*(clusterdata.largest_cluster-1)))
sys.exit(1)
# get atomic coordinates from the file (they will be shifted to the 0th cell)
array_of_lines_zero = clusters.get_0cell_coordinates(infile,clusterdata.atoms_per_molecule,clusterdata.abc)
natoms = len(array_of_lines_zero)
if (natoms%clusterdata.atoms_per_molecule != 0): "Something is wrong: natoms = %10d" % natoms
clusterdata.nmols_0cell = int(floor(natoms/clusterdata.atoms_per_molecule))
nmols_central=clusterdata.nmols_0cell
if (clusterdata.only_first_N_molecules>0):
nmols_central=clusterdata.only_first_N_molecules
clusterdata.array_of_lines = clusters.expand_0cell_to_multiple_cells(array_of_lines_zero,clusterdata.ncells,clusterdata.abc)
natoms = len(clusterdata.array_of_lines)
if (natoms%clusterdata.atoms_per_molecule != 0): "Something is wrong: natoms = %10d" % natoms
nmols = int(floor(natoms/clusterdata.atoms_per_molecule))
print "Molecules: %6d, Zero-cell molecules: %6d, Central molecules: %6d" % (nmols, clusterdata.nmols_0cell, nmols_central)
connectivityFile=infile+".connectivity-"+"%.7f" % clusterdata.Rcutoff
if ( os.path.isfile(connectivityFile) ):
clusterdata.connectivity = clusters.read_connectivity_matrix(connectivityFile,nmols)
else:
clusterdata.connectivity = clusters.create_connectivity_matrix(clusterdata.array_of_lines,nmols,clusterdata.Rcutoff)
clusters.write_connectivity_matrix(connectivityFile,clusterdata.connectivity,nmols)
# initialize all bookkeeping arrays
if (clusterdata.action==clusterdata.action_submit):
clusters.init_bookkeeping_data()
if (clusterdata.action==clusterdata.action_readenergy):
clusters.extract_DFT_energies_to_one_file(clusterdata.docluster,clusterdata.snapshotdir,clusterdata.tempdir,clusterdata.indivdir,clusterdata.energyfile)
clusters.delete_epsilon_files(clusterdata.mode_epsilon_file,clusterdata.snapshotdir,clusterdata.epsilonfile)
# get charges on atoms (to compute purely electrostatic energy)
#clusters.unpack_arc(1,clusterdata.snapshotdir,clusterdata.tempdir)
#clusterdata.array_of_charges = clusters.get_atomic_charges(nmols_unique,clusterdata.snapshotdir,clusterdata.tempdir,clusterdata.indivdir)
#clusters.delete_temp_dir(1,clusterdata.snapshotdir,clusterdata.tempdir)
# =============================================================
# main loop over clusters
# =============================================================
if (clusterdata.action==clusterdata.action_submit):
# loop over clusters of all possible sizes
cluster_index=[]
recursive.loop_over_all_clusters(cluster_index, 0, nmols_central, nmols, clusterdata.largest_cluster)
elif (clusterdata.action==clusterdata.action_readenergy or clusterdata.action == clusterdata.action_fileprep):
# loop over clusters of one particular size (select sizes by setting flags in docluster array)
for imax in range(clusterdata.largest_cluster):
if (clusterdata.docluster[imax]==1):
cluster_index=[]
strIcluster = "%03d" % (imax+1)
dbfname_Icluster = clusterdata.snapshotdir + "/" + clusterdata.dbfname + "-" + strIcluster
clusterdata.dbfhandle = open(dbfname_Icluster, 'w')
#clusterdata.dbfhandle.write("Rcutoff: ")
#clusterdata.dbfhandle.write(str(clusterdata.Rcutoff))
#clusterdata.dbfhandle.write("\n\n")
recursive.loop_over_all_clusters(cluster_index, 0, nmols_central, nmols, imax+1)
clusterdata.dbfhandle.close()
print 'Accumulated energy: %20.10f' % clusterdata.total_energy
# =============================================================
# END main loop over clusters
# =============================================================
# =============================================================
# wrap up
# =============================================================
if (clusterdata.action==clusterdata.action_submit):
clusters.close_submit_report()