6 #if .not. defined WHAM_RUN && .not. defined CLUSTER
7 use minim_data, only: maxfun,rtolf
11 #if .not. defined WHAM_RUN && .not. defined CLUSTER
12 !-----------------------------------------------------------------------------
14 !-----------------------------------------------------------------------------
15 subroutine regularize(ncart,etot,rms,cref0,iretcode)
17 use geometry, only: geom_to_var,chainbuild,var_to_geom
18 use energy, only: etotal,enerprint
19 use minimm, only: minimize
20 ! implicit real*8 (a-h,o-z)
21 ! include 'DIMENSIONS'
22 ! include 'COMMON.SBRIDGE'
23 ! include 'COMMON.CHAIN'
24 ! include 'COMMON.INTERACT'
25 ! include 'COMMON.HEADER'
26 ! include 'COMMON.IOUNITS'
27 ! include 'COMMON.MINIM'
29 real(kind=8) :: przes(3),obrot(3,3)
30 real(kind=8),dimension((nres-1)*(nres-2)/2) :: fhpb0 !(maxdim) (maxdim=(maxres-1)*(maxres-2)/2)
31 real(kind=8),dimension(6*nres) :: varia !(maxvar)(maxvar=6*maxres)
32 real(kind=8),dimension(3,ncart) :: cref0
33 real(kind=8),dimension(0:n_ene) :: energia
35 integer :: link_end0,i,maxit_reg,it
36 real(kind=8) :: etot,rms,rtolf0
37 integer :: iretcode,maxit,maxit0,maxfun0,nfun
44 print *,'Enter REGULARIZE: nnt=',nnt,' nct=',nct,' nsup=',nsup,&
45 ' nstart_seq=',nstart_seq,' nstart_sup',nstart_sup
46 write (iout,'(/a/)') 'Initial energies:'
47 call geom_to_var(nvar,varia)
48 write(iout,*) 'Warning: Calling chainbuild'
52 call enerprint(energia)
53 call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),&
54 nsup,przes,obrot,non_conv)
55 write (iout,'(a,f10.5)') &
56 'Enter REGULARIZE: Initial RMS deviation:',dsqrt(dabs(rms))
57 write (*,'(a,f10.5)') &
58 'Enter REGULARIZE: Initial RMS deviation:',dsqrt(dabs(rms))
66 print *,'Regularization: pass:',it
67 ! Minimize with distance constraints, gradually relieving the weight.
68 call minimize(etot,varia,iretcode,nfun)
70 if (iretcode.eq.11) return
71 call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),&
72 nsup,przes,obrot,non_conv)
74 write (iout,'(a,i2,a,f10.5,a,1pe14.5,a,i3/)') &
75 'Finish pass',it,', RMS deviation:',rms,', energy',etot,&
76 ' SUMSL convergence',iretcode
78 forcon(i)=0.1D0*forcon(i)
81 ! Turn off the distance constraints and re-minimize energy.
82 print *,'Final minimization ... '
86 link_end=min0(link_end,nss)
87 call minimize(etot,varia,iretcode,nfun)
89 call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),nsup,&
92 write (iout,'(a,f10.5,a,1pe14.5,a,i3/)') &
93 'Final RMS deviation:',rms,' energy',etot,' SUMSL convergence',&
99 call var_to_geom(nvar,varia)
102 end subroutine regularize
104 !-----------------------------------------------------------------------------
106 !-----------------------------------------------------------------------------
107 subroutine fitsq(rms,x,y,nn,t,b,non_conv)
109 ! implicit real*8 (a-h,o-z)
110 ! include 'COMMON.IOUNITS'
111 ! x and y are the vectors of coordinates (dimensioned (3,n)) of the two
112 ! structures to be superimposed. nn is 3*n, where n is the number of
113 ! points. t and b are respectively the translation vector and the
114 ! rotation matrix that transforms the second set of coordinates to the
115 ! frame of the first set.
116 ! eta = machine-specific variable
117 integer :: nn,i,n,j,nc
118 real(kind=8),dimension(3*nn) :: x,y
119 real(kind=8),dimension(3,3) :: b,q,r,c
120 real(kind=8),dimension(3) :: t,v,xav,yav,e
124 real(kind=8) :: rms,fn,d,sn3
127 ! small=25.0*rmdcon(3)
130 ! the following is a very lenient value for 'small'
131 real(kind=8) :: small = 0.0001D0
143 ! write(iout,*)'x = ',x(nc+i),' y = ',y(nc+i)
144 xav(i)=xav(i)+x(nc+i)/fn
145 20 yav(i)=yav(i)+y(nc+i)/fn
155 rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2
160 ! write(iout,*)'xav = ',(xav(j),j=1,3)
161 ! write(iout,*)'yav = ',(yav(j),j=1,3)
162 ! write(iout,*)'t = ',(t(j),j=1,3)
163 ! write(iout,*)'rms=',rms
164 if (rms.lt.small) return
171 rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn
173 b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn
176 call sivade(b,q,r,d,non_conv)
180 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3)
181 call mvvad(b,xav,yav,t)
184 rms=rms+2.0*c(j,i)*b(j,i)
186 if (dabs(rms).gt.small) go to 140
189 140 if (rms.gt.0.0d0) go to 150
190 ! write (iout,303) rms
193 ! 150 write (iout,302) dsqrt(rms)
196 301 format (5x,'rms deviation negligible')
197 302 format (5x,'rms deviation ',f14.6)
198 303 format (//,5x,'negative ms deviation - ',f14.6)
200 !-----------------------------------------------------------------------------
201 subroutine sivade(x,q,r,dt,non_conv)
203 ! implicit real*8(a-h,o-z)
204 ! computes q,e and r such that q(t)xr = diag(e)
205 real(kind=8),dimension(3,3) :: x,q,r
206 real(kind=8),dimension(3) :: e
207 real(kind=8),dimension(3,3) :: h,p,u
208 real(kind=8),dimension(3) :: d
212 integer :: nit,i,j,n,np,nq,npq,n0,nn
213 real(kind=8) :: dt,small,xnrm,xmax,a,den,s,b,c,dd,xn2,y,z,rt,t,v,w
216 ! write (2,*) "SIVADE"
220 ! small=2.0*rmdcon(3)
224 xnrm=xnrm+x(j,i)*x(j,i)
234 30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n))
240 den=a*(a+dabs(h(n,n)))
242 h(n,n)=h(n,n)+dsign(a,h(n,n))
249 60 x(j,i)=x(j,i)-s*h(j,n)
251 if (n.gt.1) go to 110
252 xmax=dmax1(dabs(x(1,2)),dabs(x(1,3)))
255 a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3))
256 den=a*(a+dabs(h(2,3)))
258 h(2,3)=h(2,3)+sign(a,h(2,3))
265 90 x(i,j)=x(i,j)-s*h(j,3)
270 120 p(j,i)=-d(1)*h(j,1)*h(i,1)
271 130 p(i,i)=1.0+p(i,i)
274 u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2)
275 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3)
280 ! write (2,*) "nit",nit," e",(x(i,i),i=1,3)
281 if (nit.gt.10000) then
282 print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
286 if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160
289 160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180
291 if (x(2,3).ne.0.0d0) go to 170
295 180 if (nq.eq.3) go to 310
297 ! write (2,*) "np",np," npq",npq
298 if (np.gt.npq) go to 230
302 ! write (2,*) "nn",nn
303 if (dabs(x(nn,nn)).gt.small*xnrm) go to 220
305 if (x(nn,nn+1).eq.0.0d0) go to 220
307 ! write (2,*) "nn",nn
308 go to (190,210,220),nn
310 200 call givns(x,q,1,j)
312 210 call givns(x,q,2,3)
314 ! write (2,*) "nn",nn," np",np," nq",nq," n0",n0
315 ! write (2,*) "x",(x(i,i),i=1,3)
316 if (n0.ne.0) go to 150
319 if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn)
320 b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1)
321 c=x(nn,nn)*x(nn,nn+1)
324 rt=b-xn2/(dd+sign(dsqrt(dd*dd+xn2),dd))
325 y=x(np,np)*x(np,np)-rt
326 z=x(np,np)*x(np,np+1)
328 ! write (2,*) "n",n," a",a," b",b," c",c," y",y," z",z
329 if (dabs(y).lt.dabs(z)) go to 240
331 c=1.0/dsqrt(1.0d0+t*t)
335 s=1.0/dsqrt(1.0d0+t*t)
345 260 r(j,n+1)=-s*a+c*b
348 if (dabs(y).lt.dabs(z)) go to 270
364 290 q(j,n+1)=-s*a+c*b
365 if (n.ge.nn) go to 300
375 if (nit.gt.10000) then
376 print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
380 ! write (2,*) "e",(e(i),i=1,3)
382 if (e(i).ge.0.0d0) go to 350
386 350 if (i.eq.1) go to 360
387 if (dabs(e(i)).lt.dabs(e(i-1))) go to 360
388 call switch(i,1,q,r,e)
391 if (n0.ne.0) go to 330
392 ! write (2,*) "e",(e(i),i=1,3)
393 if (dabs(e(3)).gt.small*xnrm) go to 370
395 if (dabs(e(2)).gt.small*xnrm) go to 370
397 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3))
398 ! write (2,*) "nit",nit
399 ! write (2,501) (e(i),i=1,3)
401 501 format (/,5x,'singular values - ',3e15.5)
402 end subroutine sivade
403 !-----------------------------------------------------------------------------
404 subroutine givns(a,b,m,n)
406 ! implicit real*8 (a-h,o-z)
407 real(kind=8),dimension(3,3) :: a,b
409 real(kind=8) :: t,c,s,v,w,x,y
411 if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10
430 !-----------------------------------------------------------------------------
431 subroutine switch(n,m,u,v,d)
433 ! implicit real*8 (a-h,o-z)
434 real(kind=8),dimension(3,3) :: u,v
435 real(kind=8),dimension(3) :: d
451 end subroutine switch
452 !-----------------------------------------------------------------------------
453 subroutine mvvad(b,xav,yav,t)
455 ! implicit real*8 (a-h,o-z)
456 real(kind=8),dimension(3,3) :: b
457 real(kind=8),dimension(3) :: xav,yav,t
459 ! dimension a(3,3),b(3),c(3),d(3)
463 ! 10 d(j)=d(j)+a(j,i)*b(i)
467 10 t(j)=t(j)+b(j,i)*xav(i)
470 !-----------------------------------------------------------------------------
471 real(kind=8) function det(a,b,c)
472 ! implicit real*8 (a-h,o-z)
473 real(kind=8),dimension(3) :: a,b,c
474 det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3)) &
475 +a(3)*(b(1)*c(2)-b(2)*c(1))
478 !-----------------------------------------------------------------------------
479 subroutine mmmul(a,b,c)
481 ! implicit real*8 (a-h,o-z)
482 real(kind=8),dimension(3,3) :: a,b,c
488 10 c(i,j)=c(i,j)+a(i,k)*b(k,j)
491 #if .not. defined WHAM_RUN || .not. defined CLUSTER
492 !-----------------------------------------------------------------------------
493 subroutine matvec(uvec,tmat,pvec,nback)
495 ! implicit real*8 (a-h,o-z)
496 real(kind=8),dimension(3,3) :: tmat
497 real(kind=8),dimension(3,nback) :: uvec,pvec
498 integer :: nback,i,j,k
504 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j)
507 end subroutine matvec
508 !-----------------------------------------------------------------------------
510 !-----------------------------------------------------------------------------
511 end module regularize_