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