-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhbonds.f90
179 lines (162 loc) · 4.48 KB
/
hbonds.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
! hydrogen bonding analysis for a molecule
subroutine hbonds(molx)
use parm
use atomdata, only: el
use logic, only:thresh_hbr,thresh_hba
implicit none
real(8) ang,anggrad
integer acc(6),bond(nat,nat),n_hb
type(molecule) molx
real(8) rthr,athr,a1,a2,rYX,rYH
real(8) rab,r
real(8) dummy(max_dummy)
print*,''
print*,'********************'
print*,'* H-BOND ANALYSIS *'
print*,'********************'
print*,''
print*,'type: Y-H---X'
print*,' X/Y=N,O,F,P,S,Cl'
print*,''
! bonding matrix
call bondmatrix(molx%xyz,molx%iat,bond)
! distance and angle criteria (simple!)
rthr=thresh_hbr ! Ang
athr=thresh_hba ! deg D-H--A
a1=180d0-athr
a2=180d0+athr
print*,'thresholds:'
print'(a,F8.1)', ' max Hbond length [A] :',rthr
print'(a,2(F8.1,x))',' Hbond angle range [deg] :', a1,a2
print*,''
n_hb=0
! list of active sites: N,O,F,P,S,Cl
acc(1)=7
acc(2)=8
acc(3)=9
acc(4)=15
acc(5)=16
acc(6)=17
i=0;j=0;k=0
write(*,'(a)'),' Y - H ....... X H--X [A] Y-H--X [deg] | Y-H [A] Y---X [A]'
do i=1,nat ! Y
if(.not.any(acc==molx%iat(i))) cycle
do j=1,nat ! H
if(i==j) cycle
if(bond(i,j)/=1) cycle
if(molx%iat(j)/=1) cycle
do k=1,nat ! X
if(k==i.or.k==j) cycle
if(.not.any(acc==molx%iat(k))) cycle
r=rab(molx%xyz(:,j),molx%xyz(:,k))
rYX=rab(molx%xyz(:,i),molx%xyz(:,k))
rYH=rab(molx%xyz(:,i),molx%xyz(:,j))
call angle(molx%xyz,i,j,k,ang,anggrad)
if(r<=rthr.and.anggrad>=a1.and.anggrad<=a2) then
n_hb=n_hb+1
write(*,'(a,3(x,I5,''['',a2,'']''),x,F8.4,x,F8.1,9x,2(F8.4,x))') 'H-bond', &
i,el(molx%iat(i)),&
j,el(molx%iat(j)),&
k,el(molx%iat(k)),&
r,anggrad,rYH,rYX
endif
enddo
enddo
enddo
print*,''
print'(a,I4)',' # H-bonds: ',n_hb
print*,''
end subroutine
! H-bond difference
subroutine delta_hbonds(mol1,mol2)
use parm
use atomdata, only: el
use logic, only:thresh_hbr,thresh_hba
implicit none
real(8) ang1,anggrad1
real(8) ang2,anggrad2
integer acc(6),bond(nat,nat),n_hb
type(molecule) mol1,mol2
real(8) rthr,athr,a1,a2
real(8) rab,r1,r2,da,dr,rYX,rYH,rr
real(8) dummy(max_dummy)
real(8), allocatable :: tvec(:)
print*,''
print*,'********************'
print*,'* H-BOND ANALYSIS *'
print*,'********************'
print*,''
print*,'type: Y-H---X'
print*,' X/Y=N,O,F,P,S,Cl'
print*,''
! bonding matrix
call bondmatrix(mol1%xyz,mol1%iat,bond)
! distance and angle criteria (simple!)
rthr=3.0 ! Ang
athr=120 ! deg D-H--A
rthr=thresh_hbr
athr=thresh_hba
a1=180d0-athr
a2=180d0+athr
print*,'thresholds:'
print'(a,F8.1)', ' max Hbond length [A] :',rthr
print'(a,2(F8.1,x))',' Hbond angle range [deg] :', a1,a2
print*,''
n_hb=0
! list of active sites: N,O,F,P,S,Cl
acc(1)=7
acc(2)=8
acc(3)=9
acc(4)=15
acc(5)=16
acc(6)=17
i=0;j=0;k=0
!write(*,'(a)'),' atom atom d(H--X) [A] d(Y-H--X) [deg]'
write(*,'(a)'),' Y - H ....... X H--X [A] Y-H--X [deg] | Y-H [A] Y---X [A]'
do i=1,nat ! Y
if(.not.any(acc==mol1%iat(i))) cycle
do j=1,nat ! H
if(i==j) cycle
if(bond(i,j)/=1) cycle
if(mol1%iat(j)/=1) cycle
do k=1,nat ! X
if(k==i.or.k==j) cycle
if(.not.any(acc==mol1%iat(k))) cycle
r1=rab(mol1%xyz(:,j),mol1%xyz(:,k))
rr=r1
r2=rab(mol2%xyz(:,j),mol2%xyz(:,k))
dr=r1-r2
r1=rab(mol1%xyz(:,i),mol1%xyz(:,k))
r2=rab(mol2%xyz(:,i),mol2%xyz(:,k))
rYX=r1-r2
r1=rab(mol1%xyz(:,i),mol1%xyz(:,j))
r2=rab(mol2%xyz(:,i),mol2%xyz(:,j))
rYH=r1-r2
call angle(mol1%xyz,i,j,k,ang1,anggrad1)
call angle(mol2%xyz,i,j,k,ang2,anggrad2)
da=anggrad1-anggrad2
if(rr<=rthr.and.anggrad1>=a1.and.anggrad1<=a2) then
n_hb=n_hb+1
dummy(n_hb)=dr
! write(*,'(a,x,I5,''['',a2,'']'',I5,''['',a2,'']'',x,F8.4,x,F8.1)'),'delta_H-bond',j,el(mol1%iat(i)),k,el(mol1%iat(j)),dr,da
write(*,'(a,3(x,I5,''['',a2,'']''),x,F8.4,x,F8.1,9x,2(F8.4,x))') 'H-bond', &
i,el(mol1%iat(i)),&
j,el(mol1%iat(j)),&
k,el(mol1%iat(k)),&
dr,da,rYH,rYX
endif
enddo
enddo
enddo
print*,''
print'(a,I4)',' # H-bonds: ',n_hb
print*,''
if(n_hb>0) then
allocate(tvec(n_hb))
tvec(1:n_hb)=dummy(1:n_hb)
print'(a,F8.1)',' mean abs Hbond deviation [A] : ', sum(abs(tvec))/dble(n_hb)
print'(a,F8.1)',' mean Hbond deviation [A] : ', sum(tvec)/dble(n_hb)
print'(a,F8.1)',' max (+) deviation [A] : ', maxval(tvec)
print'(a,F8.1)',' max (-) deviation [A] : ', minval(tvec)
endif
end subroutine