-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmod_input.f90
290 lines (206 loc) · 7.68 KB
/
mod_input.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
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
module mod_input
!
! Create by Veselin Kolev <vesso.kolev@gmail.com>
! 20151008233108
!
use mod_t
use mod_array
use mod_dihedral
use mod_energy
contains
subroutine readInput(files,atoms,bonds,pdihs,eps_M,sigma_M)
type(file_t), intent(in) :: files
type(atom_t), dimension(:), allocatable, intent(out) :: atoms
type(bond_t), dimension(:), allocatable, intent(out) :: bonds
type(pdih_t), dimension(:), allocatable, intent(out) :: pdihs
real(kind=real32), dimension(:,:), allocatable :: eps_M
real(kind=real32), dimension(:,:), allocatable :: sigma_M
call readAtoms(files,atoms)
atoms_s=size(atoms,1)
call readBonds(files,bonds)
bonds_s=size(bonds,1)
call readPropDihedrals(files,pdihs)
pdihs_s=size(pdihs,1)
call updateAtomAddBondsArray(atoms,bonds)
call calculateDihedralAngles(atoms,pdihs)
call calculateTorsionEnergies(pdihs)
call putAtomsInEpsGroups(atoms,eps_M)
call putAtomsInSigmaGroups(atoms,sigma_M)
end subroutine readInput
subroutine readAtoms(files,atoms)
type(file_t), intent(in) :: files
type(atom_t), dimension(:), allocatable, intent(out) :: atoms
integer(kind=int64) :: i,numAtoms
open(unit=1,file=files%atoms,form="FORMATTED")
read(1,fmt="(i10)") numAtoms
allocate(atoms(numAtoms))
do i=1,numAtoms
read(1,fmt="(i10,1x,a5,1x,i10,1x,a5,1x,a5,1x,i10,1x,f12.6,1x,&
f12.6,1x,f12.6,1x,f12.6,1x,f12.6,1x,f12.6,1x,f12.6)")&
atoms(i)%anum,atoms(i)%atype,atoms(i)%resnum,atoms(i)%resname,&
atoms(i)%aname,atoms(i)%cgrp,atoms(i)%charge,atoms(i)%mass,&
atoms(i)%sigma,atoms(i)%eps,atoms(i)%x,atoms(i)%y,atoms(i)%z
end do
close(1)
end subroutine readAtoms
subroutine readBonds(files,bonds)
type(file_t), intent(in) :: files
type(bond_t), dimension(:), allocatable, intent(out) :: bonds
integer(kind=int64) :: i,numBonds
open(unit=1,file=files%bonds,form="FORMATTED")
read(1,fmt="(i15)") numBonds
allocate(bonds(numBonds))
do i=1,numBonds
read(1,fmt="(i10,1x,i10,1x,a5,1x,a5)")&
bonds(i)%left_a,bonds(i)%right_a,&
bonds(i)%left_a_t,bonds(i)%right_a_t
end do
close(1)
end subroutine readBonds
subroutine readPropDihedrals(files,pdihs)
type(file_t), intent(in) :: files
type(pdih_t), dimension(:), allocatable, intent(out) :: pdihs
integer(kind=int64) :: i,j,numPDihrs,dummy
integer(kind=int16) :: numDihs
open(unit=1,file=files%pdihsDefs,form="FORMATTED")
open(unit=2,file=files%pdihsParams,form="FORMATTED")
read(1,fmt="(i10)") numPDihrs
allocate(pdihs(numPDihrs))
do i=1,numPDihrs
read(1,fmt="(i10,1x,i10,1x,i10,1x,i10,1x,i10,1x,i2)") &
pdihs(i)%id,&
pdihs(i)%member_id(1),pdihs(i)%member_id(2),pdihs(i)%member_id(3),&
pdihs(i)%member_id(4),numDihs
allocate(pdihs(i)%params(numDihs))
do j=1,numDihs
read(2,fmt="(i10,1x,f12.6,1x,f12.6,1x,i2)") dummy,&
pdihs(i)%params(j)%phi0,pdihs(i)%params(j)%kd,&
pdihs(i)%params(j)%mult
end do
end do
close(2)
close(1)
end subroutine readPropDihedrals
subroutine updateAtomAddBondsArray(atoms,bonds)
type(atom_t), dimension(:), intent(inout) :: atoms
type(bond_t), dimension(:), intent(in) :: bonds
integer(kind=int64) :: i
logical, dimension(size(atoms,1)) :: flag
flag(:)=.True.
do i=1,size(atoms,1)
allocate(atoms(i)%bonds(1))
end do
do i=1,size(bonds,1)
if (flag(bonds(i)%left_a)) then
flag(bonds(i)%left_a)=.False.
atoms(bonds(i)%left_a)%bonds(1)=bonds(i)%right_a
else
call extend1DIntArray(atoms(bonds(i)%left_a)%bonds,&
bonds(i)%right_a)
end if
if (flag(bonds(i)%right_a)) then
flag(bonds(i)%right_a)=.False.
atoms(bonds(i)%right_a)%bonds(1)=bonds(i)%left_a
else
call extend1DIntArray(atoms(bonds(i)%right_a)%bonds,&
bonds(i)%left_a)
end if
end do
end subroutine updateAtomAddBondsArray
subroutine putAtomsInEpsGroups(atoms,eps_M)
type(atom_t), dimension(:), intent(inout) :: atoms
real(kind=real32), dimension(:,:), allocatable, intent(out) :: eps_M
real(kind=real32), dimension(:), allocatable :: u_eps
integer(kind=int64), dimension(1) :: atoms_s
integer(kind=int64) :: i
atoms_s=shape(atoms)
allocate(u_eps(atoms_s(1)))
do i=1,atoms_s(1)
u_eps(i)=atoms(i)%eps
end do
call findUniqueElements1Darray(u_eps)
do i=1,atoms_s(1)
atoms(i)%eps_g=searchIn1DFloatArray(u_eps,atoms(i)%eps)
end do
call constructEpsMatrix(u_eps,eps_M)
end subroutine putAtomsInEpsGroups
subroutine putAtomsInSigmaGroups(atoms,sigma_M)
type(atom_t), dimension(:), intent(inout) :: atoms
real(kind=real32), dimension(:,:), allocatable, intent(out) :: &
sigma_M
real(kind=real32), dimension(:), allocatable :: u_sgm
integer(kind=int64), dimension(1) :: atoms_s
integer(kind=int64) :: i
atoms_s=shape(atoms)
allocate(u_sgm(atoms_s(1)))
do i=1,atoms_s(1)
u_sgm(i)=atoms(i)%sigma
end do
call findUniqueElements1Darray(u_sgm)
do i=1,atoms_s(1)
atoms(i)%sigma_g=searchIn1DFloatArray(u_sgm,atoms(i)%sigma)
end do
call constructSigmaMatrix(u_sgm,sigma_M)
end subroutine putAtomsInSigmaGroups
subroutine constructEpsMatrix(u_eps,eps_M)
real(kind=real32), dimension(:), intent(in) :: u_eps
real(kind=real32), dimension(:,:), allocatable, intent(out) :: eps_M
integer(kind=int64), dimension(1) :: u_eps_s
integer(kind=int64) :: i,j
u_eps_s=shape(u_eps)
allocate(eps_M(u_eps_s(1),u_eps_s(1)))
! Remarks: Do multiply the epsilon matrix elements by 4.
! It is needed by the formula used for computing VdW potential:
! V(i,j)=4.0*epsilon(i,j)*((sigma(i,j)/r(i,j))**12-(sigma(i,j)/r(i,j))**6)
do i=1,u_eps_s(1)
eps_M(i,i)=4.0_real32*u_eps(i)
do j=i+1,u_eps_s(1)
eps_M(i,j)=4.0_real32*sqrt(u_eps(i)*u_eps(j))
eps_M(j,i)=eps_M(i,j)
end do
end do
end subroutine constructEpsMatrix
subroutine constructSigmaMatrix(u_sgm,sigma_M)
real(kind=real32), dimension(:), intent(in) :: u_sgm
real(kind=real32), dimension(:,:), allocatable, intent(out) :: &
sigma_M
integer(kind=int64), dimension(1) :: u_sgm_s
integer(kind=int64) :: i,j
u_sgm_s=shape(u_sgm)
allocate(sigma_M(u_sgm_s(1),u_sgm_s(1)))
do i=1,u_sgm_s(1)
sigma_M(i,i)=u_sgm(i)
do j=i+1,u_sgm_s(1)
sigma_M(i,j)=0.5_real32*(u_sgm(i)+u_sgm(j))
sigma_M(j,i)=sigma_M(i,j)
end do
end do
end subroutine constructSigmaMatrix
subroutine createDistanceMatrix(atoms,distM)
!
! VERY IMPORTANT! To save memory we don't use two dimensional distance
! matrix. Instead an one-dimensional representation of the upper triangular
! matrix (excluding its diagonal elements) is implemented. The connection
! between the index variable of the 1D-array (q) and the indexes of the
! matrix
! elements (i and j) is:
!
! q=N*(i-1)-i*(i+1)/2+j
!
! where N is the size of the distance matrix (the number of its rows, or
! columns, or diagonal elements) - equal to the total number of the atoms.
!
type(atom_t), dimension(:), intent(in) :: atoms
real(kind=real32), dimension(:), allocatable, intent(out) :: distM
integer(kind=int64) :: i,j,counter
allocate(distM(atoms_s*(atoms_s-1)/2))
counter=1
do i=1,atoms_s-1
do j=i+1,atoms_s
distM(counter)=norm02((/atoms(i)%x,atoms(i)%y,atoms(i)%z/)-&
(/atoms(j)%x,atoms(j)%y,atoms(j)%z/))
counter=counter+1
end do
end do
end subroutine createDistanceMatrix
end module mod_input