source: metropolis.f@ bd2278d

Last change on this file since bd2278d was bd2278d, checked in by baerbaer <baerbaer@…>, 16 years ago

Reformatting comments and continuation marks.

Fortran 90 and higher use ! to mark comments no matter where they are in the
code. The only valid continuation marker is &.
I also added the SMMP.kdevelop.filelist to the repository to make it easier
to use kdevelop.

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