-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathm_fluctuations.f90
145 lines (116 loc) · 2.67 KB
/
m_fluctuations.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
!----------------------------------------------------------------------
! MELQUIADES: Metropolis Monte Carlo Program
!----------------------------------------------------------------------
!bop
!
! !Module: m_averages
!
! !Description: This module contains a routine for averages calculation.
!\\
!\\
! !Interface:
!
module m_fluctuations
!
! !Uses:
!
use m_kind
use m_simtype
use m_boxtype
use m_constants
use m_unit
use m_zeros
implicit none
real(rkind), public :: sgh, sgrho, sgvol, sget
!
! !Public members functions:
!
private
public :: r_fluctuations
!
! !Revisionn history:
! 16Nov 2015 Asdrubal Lozada
!
!eop
!---------------------------------------------
contains
!
!bop
!
! !Iroutine: r_flutuations
!
! !Description: This routine calcluates the statistical error
!
!\\
!\\
! !Interface:
!
subroutine r_fluctuations( t, y, x)
!
implicit none
!
! !Input parameters:
type(simulation), intent(inout) :: y
type(temporary), intent(inout) :: t
type(box) :: x
!
! !Revision history:
! 16Nov 2015 Asdrubal Lozada
!
!eop
!------------------------------------------------------
! Local variables
integer :: idone, k
real(rkind) :: am, cv
real(rkind) :: cnorm
open(influt,file='error',status='unknown',form='formatted')
rewind(influt)
if( .not.y%m_averag ) then
idone = 0
else
read(influt,*) idone
do k = 1, idone
read(influt,*) x%m_sehpy(k), x%m_serho(k), x%m_sevol(k), x%m_sete(k)
end do
end if
idone = idone + 1
am = 1.0_rkind / t%m_amv
!----
x%m_sehpy(idone) = (t%m_aehs*am-(t%m_aent*am)**2)/(y%m_mxmol**2)
sgh = 0.0_rkind
cv = y%m_weight * d_con / y%m_mxmol
cv = cv * cv
x%m_serho(idone) = (t%m_arhs*am -(t%m_arho*am)**2) * cv
sgrho = 0.0_rkind
cv = y%m_mxmol * d_con
cv = cv * cv
x%m_sevol(idone) = (t%m_avos*am -(t%m_avol*am)**2) / cv
sgvol = 0.0_rkind
x%m_sete(idone) = (t%m_aens*am -(t%m_aeng*am)**2)
sget = 0.0_rkind
if( y%m_averag ) then
do k = 1, idone
sgrho = sgrho + x%m_serho(k)
sget = sget + x%m_sete(k)
sgvol = sgvol + x%m_sevol(k)
sgh = sgh + x%m_sehpy(k)
end do
cnorm = real( (idone * (idone-1)), rkind )
cnorm = 1.0_rkind / cnorm
sgrho = sgrho * cnorm
sgrho = 1.0_rkind * dsqrt(dabs(sgrho))
sget = sget * cnorm
sget = 1.0_rkind * dsqrt(dabs(sget))
sgvol = sgvol * cnorm
sgvol = 1.0_rkind * dsqrt(dabs(sgvol))
sgh = sgh * cnorm
sgh = 1.0_rkind * dsqrt(dabs(sgh))
end if
rewind(influt)
write(influt,*) idone
do k = 1, idone
write(influt,*) x%m_sehpy(k), x%m_serho(k), x%m_sevol(k), x%m_sete(k)
end do
close(influt)
end subroutine r_fluctuations
end module m_fluctuations