-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmulti_bond.f90
183 lines (149 loc) · 5.02 KB
/
multi_bond.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
subroutine multi_bond
use input_parameter, &
only: jobtype, Nbeads, Nbond, TNstep, atom, atom_multi, hist_max, &
hist_min, Nfile, atom_num, Lfolding, save_beads, FNameBinary1
use calc_parameter, &
only: data_beads, data_step ! step(TNstep), beads(Nbeads,TNstep)
use calc_histogram1D, only: calc_1Dhist
implicit none
integer :: i, j, k, step
character(len=128) out_bond, out_hist
real(8) :: hist_min_ex, hist_max_ex
hist_min_ex = hist_min(1)
hist_max_ex = hist_max(1)
select case(jobtype)
case(11)
call multi_bond_all
case(12)
call multi_bond_sort
case(13:14)
call multi_bond_diff
case default
stop 'ERROR!!! wrong "Job type" option'
end select
block
integer :: Ounit
if ( save_beads .eqv. .True. ) then
open(newunit=Ounit,file=FNameBinary1, form='unformatted', access='stream', status='replace')
do step = 1, TNstep
do i = 1, Nbeads
write(Ounit) data_beads(i,step)
end do
end do
close(Ounit)
end if
end block
return
contains
! +++ jobtype == 13 +++
subroutine multi_bond_diff
integer :: Ifile
real(8) :: data_multi(Nbeads,TNstep,2), data_step_multi(TNstep,2)
write(*,*) "*****START calculating the data*****"
step = 0
do Ifile = 1, Nfile
call calc_bond_sub(Ifile,atom_num(1,Ifile),atom_num(2,Ifile),step)
end do
data_multi(:,:,1) = data_beads(:,:)
data_step_multi(:,1) = data_step(:)
step = 0
do Ifile = 1, Nfile
call calc_bond_sub(Ifile,atom_num(3,Ifile),atom_num(4,Ifile),step)
end do
data_multi(:,:,2) = data_beads(:,:)
data_step_multi(:,2) = data_step(:)
select case(jobtype)
case(13)
data_beads(:,:) = data_multi(:,:,1) - data_multi(:,:,2)
data_step(:) = data_step_multi(:,1) - data_step_multi(:,2)
write(out_hist,'("hist_diff.out")')
case(14)
data_beads(:,:) = data_multi(:,:,1) + data_multi(:,:,2)
data_step(:) = data_step_multi(:,1) + data_step_multi(:,2)
write(out_hist,'("hist_sum.out")')
end select
if ( Lfolding .eqv. .True. ) then
data_beads(:,:) = abs( data_beads(:,:) )
end if
call calc_1Dhist(out_hist_ex=out_hist)
write(out_bond, '("bond_multi_sort.out")')
open(22,file=trim(out_bond),status='replace')
do i = 1, TNstep
write(22,'(I7,10F10.5)') i, data_step(i)
end do
close(22)
end subroutine multi_bond_diff
! +++ jobtype == 11 +++
subroutine multi_bond_all
use utility
real(8) :: data_max, data_min, data_ave, data_dev
do i = 1, Nbond
data_beads(:,:) = 0.0d0
data_step(:) = 0.0d0
write(out_hist, '("hist_",a,I0,"-",a,I0".out")') &
& trim(atom(atom_multi(1,i))), atom_multi(1,i), trim(atom(atom_multi(2,i))), atom_multi(2,i)
write(out_bond, '("bond_",a,I0,"-",a,I0".out")') &
& trim(atom(atom_multi(1,i))), atom_multi(1,i), trim(atom(atom_multi(2,i))), atom_multi(2,i)
step=0
call calc_bond_sub(1,atom_multi(1,i),atom_multi(2,i),step)
data_max = maxval(data_beads)
data_min = minval(data_beads)
data_ave = sum(data_beads)/size(data_beads)
call calc_deviation(data_dev)
open(22, file=out_bond, status='replace')
write(22,*) "# ", trim(out_bond)
write(22,*) "# Maximum bond = ", data_max
write(22,*) "# Minimum bond = ", data_min
write(22,*) "# Average bond = ", data_ave
write(22,*) "# Standard deviation = ", data_dev
do k = 1, TNstep
if (mod(k,10) == 0) then
write(22,'(I7,F10.5)') k, data_step(k)
end if
! write(22,*) k, data_step(k)
end do
close(22)
call calc_1Dhist(0.0d0,0.0d0,out_hist_ex=out_hist) ! you need "data_beads"
end do
end subroutine multi_bond_all
! +++ End jobtype == 11 +++
! +++ jobtype == 12 +++
subroutine multi_bond_sort
real(8) :: data_multi(Nbeads,TNstep,Nbond), data_step_multi(TNstep,Nbond)
real(8) :: temp, temp_beads(Nbeads)
write(*,*) "*****START calculating the data*****"
do i = 1, Nbond
step = 0
call calc_bond_sub(1,atom_multi(1,i),atom_multi(2,i),step)
data_multi(:,:,i) = data_beads(:,:)
data_step_multi(:,i) = data_step(:)
end do
write(*,*) "*****START sorting*****"
do step = 1, TNstep
do i = 1, Nbond-1
do j = i+1, Nbond
if ( data_step_multi(step,i) > data_step_multi(step,j)) then
temp = data_step_multi(step,i)
data_step_multi(step,i) = data_step_multi(step,j)
data_step_multi(step,j) = temp
temp_beads(:) = data_multi(:,step,i)
data_multi(:,step,i) = data_multi(:,step,j)
data_multi(:,step,j) = temp_beads(:)
end if
end do
end do
end do
write(out_bond, '("bond_multi_sort.out")')
open(22,file=trim(out_bond),status='replace')
do i = 1, TNstep
write(22,'(I7,10F10.5)') i, data_step_multi(i,:)
end do
close(22)
do i = 1, Nbond
write(out_hist, '("hist_multisort_",I0,".out")') i
data_beads(:,:) = data_multi(:,:,i)
call calc_1Dhist(hist_min_ex, hist_max_ex, out_hist_ex=out_hist)
end do
end subroutine multi_bond_sort
! +++ End jobtype == 12 +++
end subroutine multi_bond