source/unres/regularize.F90

   1       module regularize_
   2       use names
   3       use io_units
   4       use geometry_data
   5       use energy_data
   6 #if .not. defined WHAM_RUN && .not. defined CLUSTER
   7       use minim_data, only: maxfun,rtolf
   8 #endif
   9       implicit none
  10       contains
  11 #if .not. defined WHAM_RUN && .not. defined CLUSTER
  12 !-----------------------------------------------------------------------------
  13 ! regularize.F
  14 !-----------------------------------------------------------------------------
  15       subroutine regularize(ncart,etot,rms,cref0,iretcode)
  16
  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'
  28       integer :: ncart
  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
  34       logical :: non_conv
  35       integer :: link_end0,i,maxit_reg,it
  36       real(kind=8) :: etot,rms,rtolf0
  37       integer :: iretcode,maxit,maxit0,maxfun0,nfun
  38
  39       link_end0=link_end
  40       do i=1,nhpb
  41         fhpb0(i)=forcon(i)
  42       enddo
  43       maxit_reg=2
  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'
  49       call chainbuild
  50       call etotal(energia)
  51       etot=energia(0)
  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))
  59       maxit0=maxit
  60       maxfun0=maxfun
  61       rtolf0=rtolf
  62       maxit=100
  63       maxfun=200
  64       rtolf=1.0D-2
  65       do it=1,maxit_reg
  66         print *,'Regularization: pass:',it
  67 ! Minimize with distance constraints, gradually relieving the weight.
  68         call minimize(etot,varia,iretcode,nfun)
  69         print *,'Etot=',Etot
  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)
  73         rms=dsqrt(rms)
  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
  77         do i=nss+1,nhpb
  78           forcon(i)=0.1D0*forcon(i)
  79         enddo
  80       enddo
  81 ! Turn off the distance constraints and re-minimize energy.
  82       print *,'Final minimization ... '
  83       maxit=maxit0
  84       maxfun=maxfun0
  85       rtolf=rtolf0
  86       link_end=min0(link_end,nss)
  87       call minimize(etot,varia,iretcode,nfun)
  88       print *,'Etot=',Etot
  89       call fitsq(rms,c(1,nstart_seq),cref0(1,nstart_sup-1),nsup,&
  90         przes,obrot,non_conv)
  91       rms=dsqrt(rms)
  92       write (iout,'(a,f10.5,a,1pe14.5,a,i3/)') &
  93        'Final RMS deviation:',rms,' energy',etot,' SUMSL convergence',&
  94        iretcode
  95       link_end=link_end0
  96       do i=nss+1,nhpb
  97         forcon(i)=fhpb0(i)
  98       enddo
  99       call var_to_geom(nvar,varia)
 100       call chainbuild
 101       return
 102       end subroutine regularize
 103 #endif
 104 !-----------------------------------------------------------------------------
 105 ! fitsq.f
 106 !-----------------------------------------------------------------------------
 107       subroutine fitsq(rms,x,y,nn,t,b,non_conv)
 108
 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
 121       logical :: non_conv
 122
 123 !el local variables
 124      real(kind=8) :: rms,fn,d,sn3
 125
 126 !      eta = z00100000
 127 !     small=25.0*rmdcon(3)
 128 !     small=25.0*eta
 129 !     small=25.0*10.e-10
 130 ! the following is a very lenient value for 'small'
 131       real(kind=8) :: small = 0.0001D0
 132       non_conv=.false.
 133       fn=nn
 134       do 10 i=1,3
 135       xav(i)=0.0D0
 136       yav(i)=0.0D0
 137       do 10 j=1,3
 138    10 b(j,i)=0.0D0
 139       nc=0
 140 !
 141       do 30 n=1,nn
 142       do 20 i=1,3
 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
 146    30 nc=nc+3
 147 !
 148       do i=1,3
 149         t(i)=yav(i)-xav(i)
 150       enddo
 151
 152       rms=0.0d0
 153       do n=1,nn
 154         do i=1,3
 155           rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2
 156         enddo
 157       enddo
 158       rms=dabs(rms/fn)
 159
 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
 165
 166
 167       nc=0
 168       rms=0.0D0
 169       do 50 n=1,nn
 170       do 40 i=1,3
 171       rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn
 172       do 40 j=1,3
 173       b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn
 174    40 c(j,i)=b(j,i)
 175    50 nc=nc+3
 176       call sivade(b,q,r,d,non_conv)
 177       sn3=dsign(1.0d0,d)
 178       do 120 i=1,3
 179       do 120 j=1,3
 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)
 182       do 130 i=1,3
 183       do 130 j=1,3
 184       rms=rms+2.0*c(j,i)*b(j,i)
 185   130 b(j,i)=-b(j,i)
 186       if (dabs(rms).gt.small) go to 140
 187 !     write (6,301)
 188       return
 189   140 if (rms.gt.0.0d0) go to 150
 190 !     write (iout,303) rms
 191       rms=0.0d0
 192 !     stop
 193 ! 150 write (iout,302) dsqrt(rms)
 194   150 continue
 195       return
 196   301 format (5x,'rms deviation negligible')
 197   302 format (5x,'rms deviation ',f14.6)
 198   303 format (//,5x,'negative ms deviation - ',f14.6)
 199       end subroutine fitsq
 200 !-----------------------------------------------------------------------------
 201       subroutine sivade(x,q,r,dt,non_conv)
 202
 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
 209       logical :: non_conv
 210
 211 !el local variables
 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
 214
 215 !      eta = z00100000
 216 !      write (2,*) "SIVADE"
 217       nit = 0
 218       small=25.0*10.d-10
 219 !     small=25.0*eta
 220 !     small=2.0*rmdcon(3)
 221       xnrm=0.0d0
 222       do 20 i=1,3
 223       do 10 j=1,3
 224       xnrm=xnrm+x(j,i)*x(j,i)
 225       u(j,i)=0.0d0
 226       r(j,i)=0.0d0
 227    10 h(j,i)=0.0d0
 228       u(i,i)=1.0
 229    20 r(i,i)=1.0
 230       xnrm=dsqrt(xnrm)
 231       do 110 n=1,2
 232       xmax=0.0d0
 233       do 30 j=n,3
 234    30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n))
 235       a=0.0d0
 236       do 40 j=n,3
 237       h(j,n)=x(j,n)/xmax
 238    40 a=a+h(j,n)*h(j,n)
 239       a=dsqrt(a)
 240       den=a*(a+dabs(h(n,n)))
 241       d(n)=1.0/den
 242       h(n,n)=h(n,n)+dsign(a,h(n,n))
 243       do 70 i=n,3
 244       s=0.0d0
 245       do 50 j=n,3
 246    50 s=s+h(j,n)*x(j,i)
 247       s=d(n)*s
 248       do 60 j=n,3
 249    60 x(j,i)=x(j,i)-s*h(j,n)
 250    70 continue
 251       if (n.gt.1) go to 110
 252       xmax=dmax1(dabs(x(1,2)),dabs(x(1,3)))
 253       h(2,3)=x(1,2)/xmax
 254       h(3,3)=x(1,3)/xmax
 255       a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3))
 256       den=a*(a+dabs(h(2,3)))
 257       d(3)=1.0/den
 258       h(2,3)=h(2,3)+sign(a,h(2,3))
 259       do 100 i=1,3
 260       s=0.0d0
 261       do 80 j=2,3
 262    80 s=s+h(j,3)*x(i,j)
 263       s=d(3)*s
 264       do 90 j=2,3
 265    90 x(i,j)=x(i,j)-s*h(j,3)
 266   100 continue
 267   110 continue
 268       do 130 i=1,3
 269       do 120 j=1,3
 270   120 p(j,i)=-d(1)*h(j,1)*h(i,1)
 271   130 p(i,i)=1.0+p(i,i)
 272       do 140 i=2,3
 273       do 140 j=2,3
 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)
 276       call mmmul(p,u,q)
 277   150 np=1
 278       nq=1
 279       nit=nit+1
 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!!!!!'
 283         non_conv=.true.
 284         return
 285       endif
 286       if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160
 287       x(2,3)=0.0d0
 288       nq=nq+1
 289   160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180
 290       x(1,2)=0.0d0
 291       if (x(2,3).ne.0.0d0) go to 170
 292       nq=nq+1
 293       go to 180
 294   170 np=np+1
 295   180 if (nq.eq.3) go to 310
 296       npq=4-np-nq
 297 !      write (2,*) "np",np," npq",npq
 298       if (np.gt.npq) go to 230
 299       n0=0
 300       do 220 n=np,npq
 301       nn=n+np-1
 302 !      write (2,*) "nn",nn
 303       if (dabs(x(nn,nn)).gt.small*xnrm) go to 220
 304       x(nn,nn)=0.0d0
 305       if (x(nn,nn+1).eq.0.0d0) go to 220
 306       n0=n0+1
 307 !      write (2,*) "nn",nn
 308       go to (190,210,220),nn
 309   190 do 200 j=2,3
 310   200 call givns(x,q,1,j)
 311       go to 220
 312   210 call givns(x,q,2,3)
 313   220 continue
 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
 317   230 nn=3-nq
 318       a=x(nn,nn)*x(nn,nn)
 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)
 322       dd=0.5*(a-b)
 323       xn2=c*c
 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)
 327       do 300 n=np,nn
 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
 330       t=z/y
 331       c=1.0/dsqrt(1.0d0+t*t)
 332       s=c*t
 333       go to 250
 334   240 t=y/z
 335       s=1.0/dsqrt(1.0d0+t*t)
 336       c=s*t
 337   250 do 260 j=1,3
 338       v=x(j,n)
 339       w=x(j,n+1)
 340       x(j,n)=c*v+s*w
 341       x(j,n+1)=-s*v+c*w
 342       a=r(j,n)
 343       b=r(j,n+1)
 344       r(j,n)=c*a+s*b
 345   260 r(j,n+1)=-s*a+c*b
 346       y=x(n,n)
 347       z=x(n+1,n)
 348       if (dabs(y).lt.dabs(z)) go to 270
 349       t=z/y
 350       c=1.0/dsqrt(1.0+t*t)
 351       s=c*t
 352       go to 280
 353   270 t=y/z
 354       s=1.0/dsqrt(1.0+t*t)
 355       c=s*t
 356   280 do 290 j=1,3
 357       v=x(n,j)
 358       w=x(n+1,j)
 359       a=q(j,n)
 360       b=q(j,n+1)
 361       x(n,j)=c*v+s*w
 362       x(n+1,j)=-s*v+c*w
 363       q(j,n)=c*a+s*b
 364   290 q(j,n+1)=-s*a+c*b
 365       if (n.ge.nn) go to 300
 366       y=x(n,n+1)
 367       z=x(n,n+2)
 368   300 continue
 369       go to 150
 370   310 do 320 i=1,3
 371   320 e(i)=x(i,i)
 372       nit=0
 373   330 n0=0
 374       nit=nit+1
 375       if (nit.gt.10000) then
 376         print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!'
 377         non_conv=.true.
 378         return
 379       endif
 380 !      write (2,*) "e",(e(i),i=1,3)
 381       do 360 i=1,3
 382       if (e(i).ge.0.0d0) go to 350
 383       e(i)=-e(i)
 384       do 340 j=1,3
 385   340 q(j,i)=-q(j,i)
 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)
 389       n0=n0+1
 390   360 continue
 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
 394       e(3)=0.0d0
 395       if (dabs(e(2)).gt.small*xnrm) go to 370
 396       e(2)=0.0d0
 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)
 400       return
 401   501 format (/,5x,'singular values - ',3e15.5)
 402       end subroutine sivade
 403 !-----------------------------------------------------------------------------
 404       subroutine givns(a,b,m,n)
 405
 406 !      implicit real*8 (a-h,o-z)
 407       real(kind=8),dimension(3,3) :: a,b
 408       integer :: m,n,j
 409       real(kind=8) :: t,c,s,v,w,x,y
 410
 411       if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10
 412       t=a(n,n)/a(m,n)
 413       s=1.0/dsqrt(1.0+t*t)
 414       c=s*t
 415       go to 20
 416    10 t=a(m,n)/a(n,n)
 417       c=1.0/dsqrt(1.0+t*t)
 418       s=c*t
 419    20 do 30 j=1,3
 420       v=a(m,j)
 421       w=a(n,j)
 422       x=b(j,m)
 423       y=b(j,n)
 424       a(m,j)=c*v-s*w
 425       a(n,j)=s*v+c*w
 426       b(j,m)=c*x-s*y
 427    30 b(j,n)=s*x+c*y
 428       return
 429       end subroutine givns
 430 !-----------------------------------------------------------------------------
 431       subroutine switch(n,m,u,v,d)
 432
 433 !      implicit real*8 (a-h,o-z)
 434       real(kind=8),dimension(3,3) :: u,v
 435       real(kind=8),dimension(3) :: d
 436       integer :: n,m,i
 437       real(kind=8) :: tem
 438       do 10 i=1,3
 439       tem=u(i,n)
 440       u(i,n)=u(i,n-1)
 441       u(i,n-1)=tem
 442       if (m.eq.0) go to 10
 443       tem=v(i,n)
 444       v(i,n)=v(i,n-1)
 445       v(i,n-1)=tem
 446    10 continue
 447       tem=d(n)
 448       d(n)=d(n-1)
 449       d(n-1)=tem
 450       return
 451       end subroutine switch
 452 !-----------------------------------------------------------------------------
 453       subroutine mvvad(b,xav,yav,t)
 454
 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
 458       integer :: i,j
 459 !     dimension a(3,3),b(3),c(3),d(3)
 460 !     do 10 j=1,3
 461 !     d(j)=c(j)
 462 !     do 10 i=1,3
 463 !  10 d(j)=d(j)+a(j,i)*b(i)
 464       do 10 j=1,3
 465       t(j)=yav(j)
 466       do 10 i=1,3
 467    10 t(j)=t(j)+b(j,i)*xav(i)
 468       return
 469       end subroutine mvvad
 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))
 476       return
 477       end function det
 478 !-----------------------------------------------------------------------------
 479       subroutine mmmul(a,b,c)
 480
 481 !      implicit real*8 (a-h,o-z)
 482       real(kind=8),dimension(3,3) :: a,b,c
 483       integer :: i,j,k
 484       do 10 i=1,3
 485       do 10 j=1,3
 486       c(i,j)=0.0d0
 487       do 10 k=1,3
 488    10 c(i,j)=c(i,j)+a(i,k)*b(k,j)
 489       return
 490       end subroutine mmmul
 491 #if .not. defined WHAM_RUN || .not. defined CLUSTER
 492 !-----------------------------------------------------------------------------
 493       subroutine matvec(uvec,tmat,pvec,nback)
 494
 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
 499 !
 500       do 2 j=1,nback
 501          do 1 i=1,3
 502          uvec(i,j) = 0.0d0
 503          do 1 k=1,3
 504     1    uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j)
 505     2 continue
 506       return
 507       end subroutine matvec
 508 !-----------------------------------------------------------------------------
 509 #endif
 510 !-----------------------------------------------------------------------------
 511       end module regularize_