module regularize_ use names use io_units use geometry_data use energy_data #if .not. defined WHAM_RUN && .not. defined CLUSTER use minim_data, only: maxfun,rtolf #endif implicit none contains #if .not. defined WHAM_RUN && .not. defined CLUSTER !----------------------------------------------------------------------------- ! regularize.F !----------------------------------------------------------------------------- subroutine regularize(ncart,etot,rms,cref0,iretcode) use geometry, only: geom_to_var,chainbuild,var_to_geom use energy, only: etotal,enerprint use minimm, only: minimize ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.HEADER' ! include 'COMMON.IOUNITS' ! include 'COMMON.MINIM' integer :: ncart real(kind=8) :: przes(3),obrot(3,3) real(kind=8),dimension((nres-1)*(nres-2)/2) :: fhpb0 !(maxdim) (maxdim=(maxres-1)*(maxres-2)/2) real(kind=8),dimension(6*nres) :: varia !(maxvar)(maxvar=6*maxres) real(kind=8),dimension(3,ncart) :: cref0 real(kind=8),dimension(0:n_ene) :: energia logical :: non_conv integer :: link_end0,i,maxit_reg,it real(kind=8) :: etot,rms,rtolf0 integer :: iretcode,maxit,maxit0,maxfun0,nfun link_end0=link_end do i=1,nhpb fhpb0(i)=forcon(i) enddo maxit_reg=2 print *,'Enter REGULARIZE: nnt=',nnt,' nct=',nct,' nsup=',nsup,& ' nstart_seq=',nstart_seq,' nstart_sup',nstart_sup write (iout,'(/a/)') 'Initial energies:' call geom_to_var(nvar,varia) call chainbuild call etotal(energia) etot=energia(0) call enerprint(energia) call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),& nsup,przes,obrot,non_conv) write (iout,'(a,f10.5)') & 'Enter REGULARIZE: Initial RMS deviation:',dsqrt(dabs(rms)) write (*,'(a,f10.5)') & 'Enter REGULARIZE: Initial RMS deviation:',dsqrt(dabs(rms)) maxit0=maxit maxfun0=maxfun rtolf0=rtolf maxit=100 maxfun=200 rtolf=1.0D-2 do it=1,maxit_reg print *,'Regularization: pass:',it ! Minimize with distance constraints, gradually relieving the weight. call minimize(etot,varia,iretcode,nfun) print *,'Etot=',Etot if (iretcode.eq.11) return call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),& nsup,przes,obrot,non_conv) rms=dsqrt(rms) write (iout,'(a,i2,a,f10.5,a,1pe14.5,a,i3/)') & 'Finish pass',it,', RMS deviation:',rms,', energy',etot,& ' SUMSL convergence',iretcode do i=nss+1,nhpb forcon(i)=0.1D0*forcon(i) enddo enddo ! Turn off the distance constraints and re-minimize energy. print *,'Final minimization ... ' maxit=maxit0 maxfun=maxfun0 rtolf=rtolf0 link_end=min0(link_end,nss) call minimize(etot,varia,iretcode,nfun) print *,'Etot=',Etot call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),nsup,& przes,obrot,non_conv) rms=dsqrt(rms) write (iout,'(a,f10.5,a,1pe14.5,a,i3/)') & 'Final RMS deviation:',rms,' energy',etot,' SUMSL convergence',& iretcode link_end=link_end0 do i=nss+1,nhpb forcon(i)=fhpb0(i) enddo call var_to_geom(nvar,varia) call chainbuild return end subroutine regularize #endif !----------------------------------------------------------------------------- ! fitsq.f !----------------------------------------------------------------------------- subroutine fitsq(rms,x,y,nn,t,b,non_conv) ! implicit real*8 (a-h,o-z) ! include 'COMMON.IOUNITS' ! x and y are the vectors of coordinates (dimensioned (3,n)) of the two ! structures to be superimposed. nn is 3*n, where n is the number of ! points. t and b are respectively the translation vector and the ! rotation matrix that transforms the second set of coordinates to the ! frame of the first set. ! eta = machine-specific variable integer :: nn,i,n,j,nc real(kind=8),dimension(3*nn) :: x,y real(kind=8),dimension(3,3) :: b,q,r,c real(kind=8),dimension(3) :: t,v,xav,yav,e logical :: non_conv !el local variables real(kind=8) :: rms,fn,d,sn3 ! eta = z00100000 ! small=25.0*rmdcon(3) ! small=25.0*eta ! small=25.0*10.e-10 ! the following is a very lenient value for 'small' real(kind=8) :: small = 0.0001D0 non_conv=.false. fn=nn do 10 i=1,3 xav(i)=0.0D0 yav(i)=0.0D0 do 10 j=1,3 10 b(j,i)=0.0D0 nc=0 ! do 30 n=1,nn do 20 i=1,3 ! write(iout,*)'x = ',x(nc+i),' y = ',y(nc+i) xav(i)=xav(i)+x(nc+i)/fn 20 yav(i)=yav(i)+y(nc+i)/fn 30 nc=nc+3 ! do i=1,3 t(i)=yav(i)-xav(i) enddo rms=0.0d0 do n=1,nn do i=1,3 rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2 enddo enddo rms=dabs(rms/fn) ! write(iout,*)'xav = ',(xav(j),j=1,3) ! write(iout,*)'yav = ',(yav(j),j=1,3) ! write(iout,*)'t = ',(t(j),j=1,3) ! write(iout,*)'rms=',rms if (rms.lt.small) return nc=0 rms=0.0D0 do 50 n=1,nn do 40 i=1,3 rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn do 40 j=1,3 b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn 40 c(j,i)=b(j,i) 50 nc=nc+3 call sivade(b,q,r,d,non_conv) sn3=dsign(1.0d0,d) do 120 i=1,3 do 120 j=1,3 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3) call mvvad(b,xav,yav,t) do 130 i=1,3 do 130 j=1,3 rms=rms+2.0*c(j,i)*b(j,i) 130 b(j,i)=-b(j,i) if (dabs(rms).gt.small) go to 140 ! write (6,301) return 140 if (rms.gt.0.0d0) go to 150 ! write (iout,303) rms rms=0.0d0 ! stop ! 150 write (iout,302) dsqrt(rms) 150 continue return 301 format (5x,'rms deviation negligible') 302 format (5x,'rms deviation ',f14.6) 303 format (//,5x,'negative ms deviation - ',f14.6) end subroutine fitsq !----------------------------------------------------------------------------- subroutine sivade(x,q,r,dt,non_conv) ! implicit real*8(a-h,o-z) ! computes q,e and r such that q(t)xr = diag(e) real(kind=8),dimension(3,3) :: x,q,r real(kind=8),dimension(3) :: e real(kind=8),dimension(3,3) :: h,p,u real(kind=8),dimension(3) :: d logical :: non_conv !el local variables integer :: nit,i,j,n,np,nq,npq,n0,nn real(kind=8) :: dt,small,xnrm,xmax,a,den,s,b,c,dd,xn2,y,z,rt,t,v,w ! eta = z00100000 ! write (2,*) "SIVADE" nit = 0 small=25.0*10.d-10 ! small=25.0*eta ! small=2.0*rmdcon(3) xnrm=0.0d0 do 20 i=1,3 do 10 j=1,3 xnrm=xnrm+x(j,i)*x(j,i) u(j,i)=0.0d0 r(j,i)=0.0d0 10 h(j,i)=0.0d0 u(i,i)=1.0 20 r(i,i)=1.0 xnrm=dsqrt(xnrm) do 110 n=1,2 xmax=0.0d0 do 30 j=n,3 30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n)) a=0.0d0 do 40 j=n,3 h(j,n)=x(j,n)/xmax 40 a=a+h(j,n)*h(j,n) a=dsqrt(a) den=a*(a+dabs(h(n,n))) d(n)=1.0/den h(n,n)=h(n,n)+dsign(a,h(n,n)) do 70 i=n,3 s=0.0d0 do 50 j=n,3 50 s=s+h(j,n)*x(j,i) s=d(n)*s do 60 j=n,3 60 x(j,i)=x(j,i)-s*h(j,n) 70 continue if (n.gt.1) go to 110 xmax=dmax1(dabs(x(1,2)),dabs(x(1,3))) h(2,3)=x(1,2)/xmax h(3,3)=x(1,3)/xmax a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3)) den=a*(a+dabs(h(2,3))) d(3)=1.0/den h(2,3)=h(2,3)+sign(a,h(2,3)) do 100 i=1,3 s=0.0d0 do 80 j=2,3 80 s=s+h(j,3)*x(i,j) s=d(3)*s do 90 j=2,3 90 x(i,j)=x(i,j)-s*h(j,3) 100 continue 110 continue do 130 i=1,3 do 120 j=1,3 120 p(j,i)=-d(1)*h(j,1)*h(i,1) 130 p(i,i)=1.0+p(i,i) do 140 i=2,3 do 140 j=2,3 u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2) 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3) call mmmul(p,u,q) 150 np=1 nq=1 nit=nit+1 ! write (2,*) "nit",nit," e",(x(i,i),i=1,3) if (nit.gt.10000) then print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!' non_conv=.true. return endif if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160 x(2,3)=0.0d0 nq=nq+1 160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180 x(1,2)=0.0d0 if (x(2,3).ne.0.0d0) go to 170 nq=nq+1 go to 180 170 np=np+1 180 if (nq.eq.3) go to 310 npq=4-np-nq ! write (2,*) "np",np," npq",npq if (np.gt.npq) go to 230 n0=0 do 220 n=np,npq nn=n+np-1 ! write (2,*) "nn",nn if (dabs(x(nn,nn)).gt.small*xnrm) go to 220 x(nn,nn)=0.0d0 if (x(nn,nn+1).eq.0.0d0) go to 220 n0=n0+1 ! write (2,*) "nn",nn go to (190,210,220),nn 190 do 200 j=2,3 200 call givns(x,q,1,j) go to 220 210 call givns(x,q,2,3) 220 continue ! write (2,*) "nn",nn," np",np," nq",nq," n0",n0 ! write (2,*) "x",(x(i,i),i=1,3) if (n0.ne.0) go to 150 230 nn=3-nq a=x(nn,nn)*x(nn,nn) if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn) b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1) c=x(nn,nn)*x(nn,nn+1) dd=0.5*(a-b) xn2=c*c rt=b-xn2/(dd+sign(dsqrt(dd*dd+xn2),dd)) y=x(np,np)*x(np,np)-rt z=x(np,np)*x(np,np+1) do 300 n=np,nn ! write (2,*) "n",n," a",a," b",b," c",c," y",y," z",z if (dabs(y).lt.dabs(z)) go to 240 t=z/y c=1.0/dsqrt(1.0d0+t*t) s=c*t go to 250 240 t=y/z s=1.0/dsqrt(1.0d0+t*t) c=s*t 250 do 260 j=1,3 v=x(j,n) w=x(j,n+1) x(j,n)=c*v+s*w x(j,n+1)=-s*v+c*w a=r(j,n) b=r(j,n+1) r(j,n)=c*a+s*b 260 r(j,n+1)=-s*a+c*b y=x(n,n) z=x(n+1,n) if (dabs(y).lt.dabs(z)) go to 270 t=z/y c=1.0/dsqrt(1.0+t*t) s=c*t go to 280 270 t=y/z s=1.0/dsqrt(1.0+t*t) c=s*t 280 do 290 j=1,3 v=x(n,j) w=x(n+1,j) a=q(j,n) b=q(j,n+1) x(n,j)=c*v+s*w x(n+1,j)=-s*v+c*w q(j,n)=c*a+s*b 290 q(j,n+1)=-s*a+c*b if (n.ge.nn) go to 300 y=x(n,n+1) z=x(n,n+2) 300 continue go to 150 310 do 320 i=1,3 320 e(i)=x(i,i) nit=0 330 n0=0 nit=nit+1 if (nit.gt.10000) then print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!' non_conv=.true. return endif ! write (2,*) "e",(e(i),i=1,3) do 360 i=1,3 if (e(i).ge.0.0d0) go to 350 e(i)=-e(i) do 340 j=1,3 340 q(j,i)=-q(j,i) 350 if (i.eq.1) go to 360 if (dabs(e(i)).lt.dabs(e(i-1))) go to 360 call switch(i,1,q,r,e) n0=n0+1 360 continue if (n0.ne.0) go to 330 ! write (2,*) "e",(e(i),i=1,3) if (dabs(e(3)).gt.small*xnrm) go to 370 e(3)=0.0d0 if (dabs(e(2)).gt.small*xnrm) go to 370 e(2)=0.0d0 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3)) ! write (2,*) "nit",nit ! write (2,501) (e(i),i=1,3) return 501 format (/,5x,'singular values - ',3e15.5) end subroutine sivade !----------------------------------------------------------------------------- subroutine givns(a,b,m,n) ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3,3) :: a,b integer :: m,n,j real(kind=8) :: t,c,s,v,w,x,y if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10 t=a(n,n)/a(m,n) s=1.0/dsqrt(1.0+t*t) c=s*t go to 20 10 t=a(m,n)/a(n,n) c=1.0/dsqrt(1.0+t*t) s=c*t 20 do 30 j=1,3 v=a(m,j) w=a(n,j) x=b(j,m) y=b(j,n) a(m,j)=c*v-s*w a(n,j)=s*v+c*w b(j,m)=c*x-s*y 30 b(j,n)=s*x+c*y return end subroutine givns !----------------------------------------------------------------------------- subroutine switch(n,m,u,v,d) ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3,3) :: u,v real(kind=8),dimension(3) :: d integer :: n,m,i real(kind=8) :: tem do 10 i=1,3 tem=u(i,n) u(i,n)=u(i,n-1) u(i,n-1)=tem if (m.eq.0) go to 10 tem=v(i,n) v(i,n)=v(i,n-1) v(i,n-1)=tem 10 continue tem=d(n) d(n)=d(n-1) d(n-1)=tem return end subroutine switch !----------------------------------------------------------------------------- subroutine mvvad(b,xav,yav,t) ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3,3) :: b real(kind=8),dimension(3) :: xav,yav,t integer :: i,j ! dimension a(3,3),b(3),c(3),d(3) ! do 10 j=1,3 ! d(j)=c(j) ! do 10 i=1,3 ! 10 d(j)=d(j)+a(j,i)*b(i) do 10 j=1,3 t(j)=yav(j) do 10 i=1,3 10 t(j)=t(j)+b(j,i)*xav(i) return end subroutine mvvad !----------------------------------------------------------------------------- real(kind=8) function det(a,b,c) ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3) :: a,b,c det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3)) & +a(3)*(b(1)*c(2)-b(2)*c(1)) return end function det !----------------------------------------------------------------------------- subroutine mmmul(a,b,c) ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3,3) :: a,b,c integer :: i,j,k do 10 i=1,3 do 10 j=1,3 c(i,j)=0.0d0 do 10 k=1,3 10 c(i,j)=c(i,j)+a(i,k)*b(k,j) return end subroutine mmmul #if .not. defined WHAM_RUN || .not. defined CLUSTER !----------------------------------------------------------------------------- subroutine matvec(uvec,tmat,pvec,nback) ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3,3) :: tmat real(kind=8),dimension(3,nback) :: uvec,pvec integer :: nback,i,j,k ! do 2 j=1,nback do 1 i=1,3 uvec(i,j) = 0.0d0 do 1 k=1,3 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j) 2 continue return end subroutine matvec !----------------------------------------------------------------------------- #endif !----------------------------------------------------------------------------- end module regularize_