-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfunctions.f90
227 lines (189 loc) · 5.96 KB
/
functions.f90
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
module qube_functions
use qube_error
use qube_extract
implicit none
contains
!* Select the right subroutine and execute it
subroutine gete(nmol,imet,molname,e)
implicit none
character(len=120) :: molname(*)
real*8 :: e(*)
real*8 :: e1(nmol)
real*8 :: e2(nmol)
real*8 :: edum,edum1,edum2,edum3
integer :: nmol,imet
integer :: i
character(len=120) :: aa
! common /fit/ fak1,fak2
real*8 :: fak1,fak2
logical :: vdw,ex,cp,abc,zpe,etherm
vdw=.false.
cp=.false.
abc=.false. !vdw-three-body only
zpe=.false. !add crystal zpe from (..)/freq/crystal.out
etherm=.false. !add crystal ethermal from (..)/freq/crystal.out
if(imet>=2000)then
imet=imet-2000
zpe=.true.
etherm=.true.
endif
if(imet>=1000)then
imet=imet-1000
zpe=.true.
endif
if(imet>=400)then
imet=imet-400
abc=.true.
elseif(imet>=300)then
imet=imet-300
vdw=.true.
cp=.true.
elseif(imet>=200)then
imet=imet-200
cp=.true.
elseif(imet>=100)then
imet=imet-100
vdw=.true.
endif
do i=1,nmol
e(i) =0.0d0
e1(i) =0.0d0
e2(i) =0.0d0
edum =0.0d0
edum1=0.0d0
edum2=0.0d0
!* Basic extraction
if(imet==13)then
! orca
write(aa,'(a,''/orca.out'')')trim(molname(i))
call get_en_orca(e(i),aa)
! gcp
elseif(imet==30)then
write(aa,'(a,''/.CPC'')')trim(molname(i))
call get_en_d4(e(i),aa)
! d3
elseif(imet==32) then
write(aa,'(a,''/.EDISP'')')trim(molname(i))
call get_en_d4(e(i),aa)
! total energy adf
elseif(imet==52) then
write(aa,'(a,''/adf.out'')')trim(molname(i))
call get_en_adf(e(i),aa)
! total energy dft-c
elseif(imet==53) then
write(aa,'(a,''/.dftc'')')trim(molname(i))
call get_en_d4(e(i),aa)
! nl energy from orca output
elseif(imet==54) then
write(aa,'(a,''/orca.out'')')trim(molname(i))
call get_en_orca_nl(e(i),aa)
! ri-mp2 correlation energy from orca output
elseif(imet==55) then
write(aa,'(a,''/orca.out'')')trim(molname(i))
call get_en_orca_rimp2c(e(i),aa)
! dft correlation energy from orca output
elseif(imet==56) then
write(aa,'(a,''/orca.out'')')trim(molname(i))
call get_en_orca_dftc(e(i),aa)
! ccsd(t) correlation energy from orca output
elseif(imet==57) then
write(aa,'(a,''/orca.out'')')trim(molname(i))
call get_en_orca_ccsdtc(e(i),aa)
! qchem final single point
elseif(imet==58) then
write(aa,'(a,''/qchem.out'')')trim(molname(i))
call get_en_qchem(e(i),aa)
! qchem dft correlation
elseif(imet==59) then
write(aa,'(a,''/qchem.out'')')trim(molname(i))
call get_en_qchem_dftc(e(i),aa)
! qchem ri-mos mp2 correlation
elseif(imet==60) then
write(aa,'(a,''/qchem.out'')')trim(molname(i))
call get_en_qchem_mosmp2(e(i),aa)
else
call raise('e', 'Unknown method')
endif
!* Additional corrections
! if(abc)then
! write(aa,'(a,''/dftd3.out'')')trim(molname(i))
! call getenabc(edum3,aa)
! e(i)=e(i)+edum3
! endif
if(vdw)then
write(aa,'(a,''/.edisp'')')trim(molname(i))
call get_en_d4(edum,aa)
e(i)=e(i)+edum
endif
if(cp) then
write(aa,'(a,''/.cpc'')')trim(molname(i))
call get_en_d4(edum2,aa)
e(i)=e(i)+edum2
endif
enddo
end
!> This subroutine reads command line arguments and extracts molecular names and stoch. factors.
! The molecular names are stored in the input array molname and the factors are stored in the output array fac.
! If the number of molecular names is greater than the number of factors, the subroutine throws an error.
! @param nmol The number of molecular names.
! @param molname The array of molecular names.
! @param fac The array of factors.
subroutine getreak(nmol,molname,fac)
character(len=120) :: molname(*)
character(len=120),allocatable :: dir(:)
character(len=10) tmp
real*8 :: xx(20)
real*8 :: fac(*)
logical da
integer i,j,n,nn,nfac,nmol,nargs
nargs=command_argument_count()
allocate(dir(nargs))
dir = ''
n=0
do i=1,nargs
call getarg(i,dir(i))
if(dir(i)=='x'.and.n==0)then
n=i
endif
enddo
if(n==0)n=20
nmol=0
nfac=0
do i=1,nargs
if(i<n)then
nmol=nmol+1
molname(nmol)=dir(i)
else
call readl(dir(i),xx,nn)
if(nn>0) then
nfac=nfac+1
fac(nfac)=xx(1)
endif
endif
enddo
if(nmol>nfac) error stop 'nmol > nfac'
end subroutine getreak
!> Prints the usage and available methods for the tmer2 program.
! The available methods are printed with their corresponding number.
subroutine printopt
print'(a)','usage: '//'tmer2 <dirs> ... x <stfacts> ... <imethod> [refval]',''
print'(a)','list of available <imethod>:',''
print'(6xa)','52=adf energy (adf.out)',''
print'(6xa)','53=dftc (.dftc)'
print'(6xa)','30=cpc (.cpc)'
print'(6xa)','32=d3/d4 (.edisp)',''
print'(6xa)','13=orca (orca.out)'
print'(6xa)','35=sc+nl (orca.out)'
print'(6xa)','54=nl (orca.out)'
print'(6xa)','55=ri-mp2 correlation (orca.out)'
print'(6xa)','56=dft correlation (orca.out)'
print'(6xa)','57=ccsd(t) correlation (orca.out)',''
print'(6xa)','58=total energy in the final basis set (qchem.out)'
print'(6xa)','59=dft correlation energy (qchem.out)'
print'(6xa)','60=total mos-mp2 correlation energy (qchem.out)',''
print'(6xa)','imethod+100 reads dispersion energy from .edisp'
print'(6xa)','imethod+200 reads counterpoise correcture from .cpc'
print'(6xa)','imethod+300 reads .edisp and .gcp',''
call terminate(0)
end
end module qube_functions