source/unres/regularize.F90

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