source: metropolis.f

Last change on this file was 38d77eb, checked in by baerbaer <baerbaer@…>, 14 years ago

Redirected standard out to logString.

SMMP produced a lot of log messages. This became an issue when run in massively
parallel environments. I replaced all writes to standard out to a write to logString.
The next step is to pass this string to a function that writes the messages to a log
file according to their importance and the chosen log level.

git-svn-id: ​svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@34 26dc1dd8-5c4e-0410-9ffe-d298b4865968
  • Property mode set to 100644
File size: 6.4 KB
Line 
1! **************************************************************
2!
3! This file contains the subroutines: metropolis
4!
5! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
6! Shura Hayryan, Chin-Ku
7! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
8! Jan H. Meinke, Sandipan Mohanty
9!
10! **************************************************************
11
12
13 subroutine metropolis(eol1,acz,dummy)
14!
15! SUBROUTINE FOR METROPOLIS UPDATE OF CONFIGURATIONS
16!
17! CALLS: energy,addang,grnd,dummy (function provided as argument)
18!
19 include 'INCL.H'
20 include 'INCP.H'
21 include 'incl_lund.h'
22
23 double precision eol, energy, grnd, dv, addang, enw, delta, dummy
24 double precision rd, ex, acz, eol1
25
26 integer nsw, iupstate, iupt, ivar, jv, iml, i, ncalls, nacalls
27
28 external dummy
29! common/bet/beta
30 common/updstats/ncalls(5),nacalls(5)
31 integer updtch1, updtch2,bgs
32 double precision vrol(mxvr), gbol
33!f2py intent(in, out) eol
34!f2py intent(in, out) acz
35 eol = energy()
36! external rand
37 do nsw=1,nvr
38! Loop over dihedrals
39!
40 iupstate=0
41 iupt=1
42 if (rndord) then
43 ivar = 1+int(nvr*grnd())
44 else
45 ivar=nsw
46 endif
47 jv = idvr(ivar)
48 if (nmvr(jv)(1:1).eq.'x') then
49 iupt=1
50 else
51 if (upchswitch.eq.0) then
52 iupt=1
53 else if (upchswitch.eq.1) then
54 iupt=updtch1(ivar,beta)
55 else
56 iupt=updtch2(ivar,beta)
57 endif
58 endif
59 if (iupt.eq.2) then
60 iupstate=bgs(eol,dummy)
61 else
62! Simple twist of
63! Get Proposal configuration
64 vrol=vlvr!(jv)
65 dv=axvr(jv)*(grnd()-0.5)
66 vlvr(jv)=addang(vrol(jv),dv)
67!
68! Get dummy of proposal configuration
69!
70 enw = energy()
71!
72 delta = dummy(enw) - dummy(eol)
73! ___________________________ check acceptance criteria
74 if (delta.le.0.0d0) then
75 eol=enw
76 iupstate=1
77 else
78 rd=grnd()
79 delta = min(delta,100.0d0)
80 ex=exp( -delta )
81 if ( ex .ge. rd ) then
82 eol=enw
83 iupstate=1
84 else
85 vlvr=vrol
86 endif
87 endif
88 endif
89 acz=acz+iupstate*1.0
90 call accanalyze(iupt,iupstate)
91 end do
92
93! Updates on relative position of different molecules when there are many
94 if (ntlml.gt.1) then
95 do iml=1, ntlml
96 do i=1,3
97 iupstate=0
98 gbol = gbpr(i, iml)
99 gbpr(i, iml) = gbpr(i, iml) + (grnd()-0.5)
100 if (boxsize.gt.0) then
101 if (gbpr(i, iml).gt.boxsize) then
102 gbpr(i, iml)=boxsize
103 elseif (gbpr(i, iml).lt.-boxsize) then
104 gbpr(i, iml)=-boxsize
105 endif
106 endif
107!
108! Get dummy of proposal configuration
109!
110 enw = energy()
111!
112 delta = dummy(enw) - dummy(eol)
113!
114! ____________________________ check acceptance criteria
115!
116 if (delta.le.0.0d0) then
117 eol=enw
118 iupstate=1
119 else
120 rd=grnd()
121 delta = min(delta,100.0d0)
122 ex=exp( -delta )
123 if ( ex .ge. rd ) then
124 eol=enw
125 iupstate=1
126 else
127 gbpr(i, iml)=gbol
128 endif
129 endif
130 acz=acz+iupstate*1.0
131 call accanalyze(3,iupstate)
132 enddo
133 do i=4,6
134 iupstate=0
135 gbol = gbpr(i, iml)
136 if (i.eq.5) then ! global psi value must be in [-pi/2, pi/2]
137 gbpr(i, iml) = (grnd()-0.5) * pi
138 else
139 gbpr(i, iml) = (grnd()-0.5) * pi2
140 endif
141!
142! Get dummy of proposal configuration
143!
144 enw = energy()
145!
146 delta = dummy(enw) - dummy(eol)
147!
148! ____________________________ check acceptance criteria
149!
150 if (delta.le.0.0d0) then
151 eol=enw
152 iupstate=1
153 else
154 rd=grnd()
155 delta = min(delta,100.0d0)
156 ex=exp( -delta )
157 if ( ex .ge. rd ) then
158 eol=enw
159 iupstate=1
160 else
161 gbpr(i, iml)=gbol
162 endif
163 endif
164 acz=acz+iupstate*1.0
165 call accanalyze(4,iupstate)
166 enddo
167 enddo
168 endif
169!
170! Re-calculate energy
171!
172 enw = energy()
173 if(abs(eol-enw).gt.0.000001) then
174 write (logString, *) 'metropolis: eol,enw,difference:', eol,
175 & enw, eol-enw
176 if (eol.lt.100000) then
177 stop 'metropolis: eol and enw difference unacceptable'
178 endif
179 endif
180!
181 eol1 = eol
182 return
183!
184 end
185
186 subroutine accanalyze(iuptype,iupdstate)
187 integer ncalls, nacalls, iupdstate, iuptype
188 common/updstats/ncalls(5),nacalls(5)
189 ncalls(5)=ncalls(5)+1
190 nacalls(5)=nacalls(5)+iupdstate
191 ncalls(iuptype)=ncalls(iuptype)+1
192 nacalls(iuptype)=nacalls(iuptype)+iupdstate
193 end
194
195
196 integer function updtch1(iiii,bbbb)
197 logical rndord
198 integer upchswitch, iiii
199 double precision bgsprob, grnd, bbbb
200 common /updchois/rndord,upchswitch,bgsprob
201 if (grnd().lt.bgsprob) then
202 updtch1=2
203 else
204 updtch1=1
205 endif
206 return
207 end
208
209 integer function updtch2(iiii,bbbb)
210 integer iiii
211 double precision bbbb,curprob, grnd,up2bmax,up2bmin
212
213 common/updtparam/up2bmax,up2bmin
214 curprob=(bbbb-up2bmin)/(up2bmax-up2bmin)
215 if (grnd().lt.curprob) then
216 updtch2=2
217 else
218 updtch2=1
219 endif
220 end
221 block data updtchs
222 integer ncalls, nacalls
223 double precision up2bmax, up2bmin
224 common/updtparam/up2bmax,up2bmin
225 data up2bmax,up2bmin/4.0d0,0.5d0/
226 common/updstats/ncalls(5),nacalls(5)
227 data ncalls/5*0/
228 data nacalls/5*0/
229 end
Note: See TracBrowser for help on using the repository browser.