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