source/unres/src-HCD-5D/brown_step.F

   1 c-------------------------------------------------------------------------------
   2       subroutine brown_step(itime)
   3 c------------------------------------------------
   4 c  Perform a single Euler integration step of Brownian dynamics
   5 c------------------------------------------------
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8 #ifdef MPI
   9       include 'mpif.h'
  10 #endif
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.VAR'
  13       include 'COMMON.MD'
  14 #ifdef FIVEDIAG
  15        include 'COMMON.LAGRANGE.5diag'
  16 #else
  17        include 'COMMON.LAGRANGE'
  18 #endif
  19 #ifndef LANG0
  20       include 'COMMON.LANGEVIN'
  21 #else
  22 #ifdef FIVEDIAG
  23       include 'COMMON.LANGEVIN.lang0.5diag'
  24 #else
  25       include 'COMMON.LANGEVIN.lang0'
  26 #endif
  27 #endif
  28       include 'COMMON.CHAIN'
  29       include 'COMMON.DERIV'
  30       include 'COMMON.GEO'
  31       include 'COMMON.LOCAL'
  32       include 'COMMON.INTERACT'
  33       include 'COMMON.IOUNITS'
  34       include 'COMMON.NAMES'
  35       include 'COMMON.TIME1'
  36       double precision zapas(MAXRES6)
  37       integer ilen,rstcount
  38       external ilen
  39       double precision stochforcvec(MAXRES6)
  40       double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
  41      & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
  42      & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
  43      & ppvec(maxres2)
  44       common /stochcalc/ stochforcvec
  45       common /gucio/ cm
  46       integer itime
  47       logical lprn /.false./,lprn1 /.false./
  48       integer maxiter /5/
  49       double precision difftol /1.0d-5/
  50       nbond=nct-nnt
  51       do i=nnt,nct
  52         if (itype(i).ne.10) nbond=nbond+1
  53       enddo
  54 c
  55       if (lprn1) then
  56         write (iout,*) "Generalized inverse of fricmat"
  57         call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
  58       endif
  59       do i=1,dimen
  60         do j=1,nbond
  61           Bmat(i,j)=0.0d0
  62         enddo
  63       enddo
  64       ind=3
  65       ind1=0
  66       do i=nnt,nct-1
  67         ind1=ind1+1
  68         do j=1,3
  69           Bmat(ind+j,ind1)=dC_norm(j,i)
  70         enddo
  71         ind=ind+3
  72       enddo
  73       do i=nnt,nct
  74         if (itype(i).ne.10) then
  75           ind1=ind1+1
  76           do j=1,3
  77             Bmat(ind+j,ind1)=dC_norm(j,i+nres)
  78           enddo
  79           ind=ind+3
  80         endif
  81       enddo
  82       if (lprn1) then
  83         write (iout,*) "Matrix Bmat"
  84         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
  85       endif
  86       do i=1,dimen
  87         do j=1,nbond
  88           GBmat(i,j)=0.0d0
  89           do k=1,dimen
  90             GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
  91           enddo
  92         enddo
  93       enddo
  94       if (lprn1) then
  95         write (iout,*) "Matrix GBmat"
  96         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
  97       endif
  98       do i=1,nbond
  99         do j=1,nbond
 100           Cmat(i,j)=0.0d0
 101           do k=1,dimen
 102             Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
 103           enddo
 104         enddo
 105       enddo
 106       if (lprn1) then
 107         write (iout,*) "Matrix Cmat"
 108         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
 109       endif
 110       call matinvert(nbond,MAXRES2,Cmat,Cinv)
 111       if (lprn1) then
 112         write (iout,*) "Matrix Cinv"
 113         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
 114       endif
 115       do i=1,dimen
 116         do j=1,nbond
 117           Tmat(i,j)=0.0d0
 118           do k=1,nbond
 119             Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
 120           enddo
 121         enddo
 122       enddo
 123       if (lprn1) then
 124         write (iout,*) "Matrix Tmat"
 125         call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
 126       endif
 127       do i=1,dimen
 128         do j=1,dimen
 129           if (i.eq.j) then
 130             Pmat(i,j)=1.0d0
 131           else
 132             Pmat(i,j)=0.0d0
 133           endif
 134           do k=1,nbond
 135             Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
 136           enddo
 137         enddo
 138       enddo
 139       if (lprn1) then
 140         write (iout,*) "Matrix Pmat"
 141         call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
 142       endif
 143       do i=1,dimen
 144         Td(i)=0.0d0
 145         ind=0
 146         do k=nnt,nct-1
 147           ind=ind+1
 148           Td(i)=Td(i)+vbl*Tmat(i,ind)
 149         enddo
 150         do k=nnt,nct
 151           if (itype(k).ne.10) then
 152             ind=ind+1
 153             Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
 154           endif
 155         enddo
 156       enddo
 157       if (lprn1) then
 158         write (iout,*) "Vector Td"
 159         do i=1,dimen
 160           write (iout,'(i5,f10.5)') i,Td(i)
 161         enddo
 162       endif
 163       call stochastic_force(stochforcvec)
 164       if (lprn) then
 165         write (iout,*) "stochforcvec"
 166         do i=1,dimen
 167           write (iout,*) i,stochforcvec(i)
 168         enddo
 169       endif
 170       do j=1,3
 171         zapas(j)=-gcart(j,0)+stochforcvec(j)
 172         d_t_work(j)=d_t(j,0)
 173         dC_work(j)=dC_old(j,0)
 174       enddo
 175       ind=3
 176       do i=nnt,nct-1
 177         do j=1,3
 178           ind=ind+1
 179           zapas(ind)=-gcart(j,i)+stochforcvec(ind)
 180           dC_work(ind)=dC_old(j,i)
 181         enddo
 182       enddo
 183       do i=nnt,nct
 184         if (itype(i).ne.10) then
 185           do j=1,3
 186             ind=ind+1
 187             zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
 188             dC_work(ind)=dC_old(j,i+nres)
 189           enddo
 190         endif
 191       enddo
 192
 193       if (lprn) then
 194         write (iout,*) "Initial d_t_work"
 195         do i=1,dimen
 196           write (iout,*) i,d_t_work(i)
 197         enddo
 198       endif
 199
 200       do i=1,dimen
 201         d_t_work(i)=0.0d0
 202         do j=1,dimen
 203           d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
 204         enddo
 205       enddo
 206
 207       do i=1,dimen
 208         zapas(i)=Td(i)
 209         do j=1,dimen
 210           zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
 211         enddo
 212       enddo
 213       if (lprn1) then
 214         write (iout,*) "Final d_t_work and zapas"
 215         do i=1,dimen
 216           write (iout,*) i,d_t_work(i),zapas(i)
 217         enddo
 218       endif
 219
 220       do j=1,3
 221         d_t(j,0)=d_t_work(j)
 222         dc(j,0)=zapas(j)
 223         dc_work(j)=dc(j,0)
 224       enddo
 225       ind=3
 226       do i=nnt,nct-1
 227         do j=1,3
 228           d_t(j,i)=d_t_work(i)
 229           dc(j,i)=zapas(ind+j)
 230           dc_work(ind+j)=dc(j,i)
 231         enddo
 232         ind=ind+3
 233       enddo
 234       do i=nnt,nct
 235         do j=1,3
 236           d_t(j,i+nres)=d_t_work(ind+j)
 237           dc(j,i+nres)=zapas(ind+j)
 238           dc_work(ind+j)=dc(j,i+nres)
 239         enddo
 240         ind=ind+3
 241       enddo
 242       if (lprn) then
 243         call chainbuild_cart
 244         write (iout,*) "Before correction for rotational lengthening"
 245         write (iout,*) "New coordinates",
 246      &  " and differences between actual and standard bond lengths"
 247         ind=0
 248         do i=nnt,nct-1
 249           ind=ind+1
 250           xx=vbld(i+1)-vbl
 251           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 252      &        i,(dC(j,i),j=1,3),xx
 253         enddo
 254         do i=nnt,nct
 255           if (itype(i).ne.10) then
 256             ind=ind+1
 257             xx=vbld(i+nres)-vbldsc0(1,itype(i))
 258             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 259      &       i,(dC(j,i+nres),j=1,3),xx
 260           endif
 261         enddo
 262       endif
 263 c Second correction (rotational lengthening)
 264 c      do iter=1,maxiter
 265       diffmax=0.0d0
 266       ind=0
 267       do i=nnt,nct-1
 268         ind=ind+1
 269         blen2 = scalar(dc(1,i),dc(1,i))
 270         ppvec(ind)=2*vbl**2-blen2
 271         diffbond=dabs(vbl-dsqrt(blen2))
 272         if (diffbond.gt.diffmax) diffmax=diffbond
 273         if (ppvec(ind).gt.0.0d0) then
 274           ppvec(ind)=dsqrt(ppvec(ind))
 275         else
 276           ppvec(ind)=0.0d0
 277         endif
 278         if (lprn) then
 279           write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
 280         endif
 281       enddo
 282       do i=nnt,nct
 283         if (itype(i).ne.10) then
 284           ind=ind+1
 285           blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
 286           ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
 287           diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
 288           if (diffbond.gt.diffmax) diffmax=diffbond
 289           if (ppvec(ind).gt.0.0d0) then
 290             ppvec(ind)=dsqrt(ppvec(ind))
 291           else
 292             ppvec(ind)=0.0d0
 293           endif
 294           if (lprn) then
 295             write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
 296           endif
 297         endif
 298       enddo
 299       if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
 300       if (diffmax.lt.difftol) goto 10
 301       do i=1,dimen
 302         Td(i)=0.0d0
 303         do j=1,nbond
 304           Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
 305         enddo
 306       enddo
 307       do i=1,dimen
 308         zapas(i)=Td(i)
 309         do j=1,dimen
 310           zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
 311         enddo
 312       enddo
 313       do j=1,3
 314         dc(j,0)=zapas(j)
 315         dc_work(j)=zapas(j)
 316       enddo
 317       ind=3
 318       do i=nnt,nct-1
 319         do j=1,3
 320           dc(j,i)=zapas(ind+j)
 321           dc_work(ind+j)=zapas(ind+j)
 322         enddo
 323         ind=ind+3
 324       enddo
 325       do i=nnt,nct
 326         if (itype(i).ne.10) then
 327           do j=1,3
 328             dc(j,i+nres)=zapas(ind+j)
 329             dc_work(ind+j)=zapas(ind+j)
 330           enddo
 331           ind=ind+3
 332         endif
 333       enddo
 334 c   Building the chain from the newly calculated coordinates
 335       call chainbuild_cart
 336       if(ntwe.ne.0) then
 337       if (large.and. mod(itime,ntwe).eq.0) then
 338         write (iout,*) "Cartesian and internal coordinates: step 1"
 339         call cartprint
 340         call intout
 341         write (iout,'(a)') "Potential forces"
 342         do i=0,nres
 343           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
 344      &    (-gxcart(j,i),j=1,3)
 345         enddo
 346         write (iout,'(a)') "Stochastic forces"
 347         do i=0,nres
 348           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
 349      &    (stochforc(j,i+nres),j=1,3)
 350         enddo
 351         write (iout,'(a)') "Velocities"
 352         do i=0,nres
 353           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 354      &    (d_t(j,i+nres),j=1,3)
 355         enddo
 356       endif
 357       endif
 358       if (lprn) then
 359         write (iout,*) "After correction for rotational lengthening"
 360         write (iout,*) "New coordinates",
 361      &  " and differences between actual and standard bond lengths"
 362         ind=0
 363         do i=nnt,nct-1
 364           ind=ind+1
 365           xx=vbld(i+1)-vbl
 366           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 367      &        i,(dC(j,i),j=1,3),xx
 368         enddo
 369         do i=nnt,nct
 370           if (itype(i).ne.10) then
 371             ind=ind+1
 372             xx=vbld(i+nres)-vbldsc0(1,itype(i))
 373             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 374      &       i,(dC(j,i+nres),j=1,3),xx
 375           endif
 376         enddo
 377       endif
 378 c      ENDDO
 379 c      write (iout,*) "Too many attempts at correcting the bonds"
 380 c      stop
 381    10 continue
 382 #ifdef MPI
 383       tt0 =MPI_Wtime()
 384 #else
 385       tt0 = tcpu()
 386 #endif
 387 c Calculate energy and forces
 388       call zerograd
 389       call etotal(potEcomp)
 390       potE=potEcomp(0)-potEcomp(27)
 391       call cartgrad
 392       totT=totT+d_time
 393       totTafm=totT
 394 c  Calculate the kinetic and total energy and the kinetic temperature
 395       call kinetic(EK)
 396 #ifdef MPI
 397       t_enegrad=t_enegrad+MPI_Wtime()-tt0
 398 #else
 399       t_enegrad=t_enegrad+tcpu()-tt0
 400 #endif
 401       totE=EK+potE
 402       kinetic_T=2.0d0/(dimen*Rb)*EK
 403       return
 404       end
 405