source/unres/src_MD-M/stochfric.F

   1       subroutine friction_force
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.VAR'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.DERIV'
   7       include 'COMMON.GEO'
   8       include 'COMMON.LOCAL'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.MD'
  11 #ifndef LANG0
  12       include 'COMMON.LANGEVIN'
  13 #else
  14       include 'COMMON.LANGEVIN.lang0'
  15 #endif
  16       include 'COMMON.IOUNITS'
  17       double precision gamvec(MAXRES6)
  18       common /syfek/ gamvec
  19       double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
  20      & ginvfric(maxres2,maxres2)
  21       common /przechowalnia/ ginvfric
  22
  23       logical lprn /.false./, checkmode /.false./
  24
  25       do i=0,MAXRES2
  26         do j=1,3
  27           friction(j,i)=0.0d0
  28         enddo
  29       enddo
  30
  31       do j=1,3
  32         d_t_work(j)=d_t(j,0)
  33       enddo
  34       ind=3
  35       do i=nnt,nct-1
  36         do j=1,3
  37           d_t_work(ind+j)=d_t(j,i)
  38         enddo
  39         ind=ind+3
  40       enddo
  41       do i=nnt,nct
  42         if (itype(i).ne.10 .and. itype(i).ne.21) then
  43           do j=1,3
  44             d_t_work(ind+j)=d_t(j,i+nres)
  45           enddo
  46           ind=ind+3
  47         endif
  48       enddo
  49
  50       call fricmat_mult(d_t_work,fric_work)
  51
  52       if (.not.checkmode) return
  53
  54       if (lprn) then
  55         write (iout,*) "d_t_work and fric_work"
  56         do i=1,3*dimen
  57           write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
  58         enddo
  59       endif
  60       do j=1,3
  61         friction(j,0)=fric_work(j)
  62       enddo
  63       ind=3
  64       do i=nnt,nct-1
  65         do j=1,3
  66           friction(j,i)=fric_work(ind+j)
  67         enddo
  68         ind=ind+3
  69       enddo
  70       do i=nnt,nct
  71         if (itype(i).ne.10 .and. itype(i).ne.21) then
  72           do j=1,3
  73             friction(j,i+nres)=fric_work(ind+j)
  74           enddo
  75           ind=ind+3
  76         endif
  77       enddo
  78       if (lprn) then
  79         write(iout,*) "Friction backbone"
  80         do i=0,nct-1
  81           write(iout,'(i5,3e15.5,5x,3e15.5)')
  82      &     i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
  83         enddo
  84         write(iout,*) "Friction side chain"
  85         do i=nnt,nct
  86           write(iout,'(i5,3e15.5,5x,3e15.5)')
  87      &     i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
  88         enddo
  89       endif
  90       if (lprn) then
  91         do j=1,3
  92           vv(j)=d_t(j,0)
  93         enddo
  94         do i=nnt,nct
  95           do j=1,3
  96             vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
  97             vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
  98             vv(j)=vv(j)+d_t(j,i)
  99           enddo
 100         enddo
 101         write (iout,*) "vvtot backbone and sidechain"
 102         do i=nnt,nct
 103           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
 104      &     (vvtot(j,i+nres),j=1,3)
 105         enddo
 106         ind=0
 107         do i=nnt,nct-1
 108           do j=1,3
 109             v_work(ind+j)=vvtot(j,i)
 110           enddo
 111           ind=ind+3
 112         enddo
 113         do i=nnt,nct
 114           do j=1,3
 115             v_work(ind+j)=vvtot(j,i+nres)
 116           enddo
 117           ind=ind+3
 118         enddo
 119         write (iout,*) "v_work gamvec and site-based friction forces"
 120         do i=1,dimen1
 121           write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
 122      &      gamvec(i)*v_work(i)
 123         enddo
 124 c        do i=1,dimen
 125 c          fric_work1(i)=0.0d0
 126 c          do j=1,dimen1
 127 c            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
 128 c          enddo
 129 c        enddo
 130 c        write (iout,*) "fric_work and fric_work1"
 131 c        do i=1,dimen
 132 c          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
 133 c        enddo
 134         do i=1,dimen
 135           do j=1,dimen
 136             ginvfric(i,j)=0.0d0
 137             do k=1,dimen
 138               ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
 139             enddo
 140           enddo
 141         enddo
 142         write (iout,*) "ginvfric"
 143         do i=1,dimen
 144           write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
 145         enddo
 146         write (iout,*) "symmetry check"
 147         do i=1,dimen
 148           do j=1,i-1
 149             write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
 150           enddo
 151         enddo
 152       endif
 153       return
 154       end
 155 c-----------------------------------------------------
 156       subroutine stochastic_force(stochforcvec)
 157       implicit real*8 (a-h,o-z)
 158       include 'DIMENSIONS'
 159 #ifdef MPI
 160       include 'mpif.h'
 161 #endif
 162       include 'COMMON.VAR'
 163       include 'COMMON.CHAIN'
 164       include 'COMMON.DERIV'
 165       include 'COMMON.GEO'
 166       include 'COMMON.LOCAL'
 167       include 'COMMON.INTERACT'
 168       include 'COMMON.MD'
 169       include 'COMMON.TIME1'
 170 #ifndef LANG0
 171       include 'COMMON.LANGEVIN'
 172 #else
 173       include 'COMMON.LANGEVIN.lang0'
 174 #endif
 175       include 'COMMON.IOUNITS'
 176
 177       double precision x,sig,lowb,highb,
 178      & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
 179      & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
 180       logical lprn /.false./
 181       do i=0,MAXRES2
 182         do j=1,3
 183           stochforc(j,i)=0.0d0
 184         enddo
 185       enddo
 186       x=0.0d0
 187
 188       time00=MPI_Wtime()
 189 c Compute the stochastic forces acting on bodies. Store in force.
 190       do i=nnt,nct-1
 191         sig=stdforcp(i)
 192         lowb=-5*sig
 193         highb=5*sig
 194         do j=1,3
 195           force(j,i)=anorm_distr(x,sig,lowb,highb)
 196         enddo
 197       enddo
 198       do i=nnt,nct
 199         sig2=stdforcsc(i)
 200         lowb2=-5*sig2
 201         highb2=5*sig2
 202         do j=1,3
 203           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
 204         enddo
 205       enddo
 206       time_fsample=time_fsample+MPI_Wtime()-time00
 207 c Compute the stochastic forces acting on virtual-bond vectors.
 208       do j=1,3
 209         ff(j)=0.0d0
 210       enddo
 211       do i=nct-1,nnt,-1
 212         do j=1,3
 213           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
 214         enddo
 215         do j=1,3
 216           ff(j)=ff(j)+force(j,i)
 217         enddo
 218         if (itype(i+1).ne.21) then
 219           do j=1,3
 220             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
 221             ff(j)=ff(j)+force(j,i+nres+1)
 222           enddo
 223         endif
 224       enddo
 225       do j=1,3
 226         stochforc(j,0)=ff(j)+force(j,nnt+nres)
 227       enddo
 228       do i=nnt,nct
 229         if (itype(i).ne.10 .and. itype(i).ne.21) then
 230           do j=1,3
 231             stochforc(j,i+nres)=force(j,i+nres)
 232           enddo
 233         endif
 234       enddo
 235
 236       do j=1,3
 237         stochforcvec(j)=stochforc(j,0)
 238       enddo
 239       ind=3
 240       do i=nnt,nct-1
 241         do j=1,3
 242           stochforcvec(ind+j)=stochforc(j,i)
 243         enddo
 244         ind=ind+3
 245       enddo
 246       do i=nnt,nct
 247         if (itype(i).ne.10 .and. itype(i).ne.21) then
 248           do j=1,3
 249             stochforcvec(ind+j)=stochforc(j,i+nres)
 250           enddo
 251           ind=ind+3
 252         endif
 253       enddo
 254       if (lprn) then
 255         write (iout,*) "stochforcvec"
 256         do i=1,3*dimen
 257           write(iout,'(i5,e15.5)') i,stochforcvec(i)
 258         enddo
 259         write(iout,*) "Stochastic forces backbone"
 260         do i=0,nct-1
 261           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
 262         enddo
 263         write(iout,*) "Stochastic forces side chain"
 264         do i=nnt,nct
 265           write(iout,'(i5,3e15.5)')
 266      &      i,(stochforc(j,i+nres),j=1,3)
 267         enddo
 268       endif
 269
 270       if (lprn) then
 271
 272       ind=0
 273       do i=nnt,nct-1
 274         write (iout,*) i,ind
 275         do j=1,3
 276           forcvec(ind+j)=force(j,i)
 277         enddo
 278         ind=ind+3
 279       enddo
 280       do i=nnt,nct
 281         write (iout,*) i,ind
 282         do j=1,3
 283           forcvec(j+ind)=force(j,i+nres)
 284         enddo
 285         ind=ind+3
 286       enddo
 287
 288       write (iout,*) "forcvec"
 289       ind=0
 290       do i=nnt,nct-1
 291         do j=1,3
 292           write (iout,'(2i3,2f10.5)') i,j,force(j,i),
 293      &      forcvec(ind+j)
 294         enddo
 295         ind=ind+3
 296       enddo
 297       do i=nnt,nct
 298         do j=1,3
 299           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
 300      &     forcvec(ind+j)
 301         enddo
 302         ind=ind+3
 303       enddo
 304
 305       endif
 306
 307       return
 308       end
 309 c------------------------------------------------------------------
 310       subroutine setup_fricmat
 311       implicit real*8 (a-h,o-z)
 312       include 'mpif.h'
 313       include 'DIMENSIONS'
 314       include 'COMMON.VAR'
 315       include 'COMMON.CHAIN'
 316       include 'COMMON.DERIV'
 317       include 'COMMON.GEO'
 318       include 'COMMON.LOCAL'
 319       include 'COMMON.INTERACT'
 320       include 'COMMON.MD'
 321       include 'COMMON.SETUP'
 322       include 'COMMON.TIME1'
 323 c      integer licznik /0/
 324 c      save licznik
 325 #ifndef LANG0
 326       include 'COMMON.LANGEVIN'
 327 #else
 328       include 'COMMON.LANGEVIN.lang0'
 329 #endif
 330       include 'COMMON.IOUNITS'
 331       integer IERROR
 332       integer i,j,ind,ind1,m
 333       logical lprn /.false./
 334       double precision dtdi,gamvec(MAXRES2),
 335      &  ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2)
 336       common /syfek/ gamvec
 337       double precision work(8*maxres2)
 338       integer iwork(maxres2)
 339       common /przechowalnia/ ginvfric,Ghalf,fcopy
 340 #ifdef MPI
 341       if (fg_rank.ne.king) goto 10
 342 #endif
 343 c  Zeroing out fricmat
 344       do i=1,dimen
 345         do j=1,dimen
 346           fricmat(i,j)=0.0d0
 347         enddo
 348       enddo
 349 c  Load the friction coefficients corresponding to peptide groups
 350       ind1=0
 351       do i=nnt,nct-1
 352         ind1=ind1+1
 353         gamvec(ind1)=gamp
 354       enddo
 355 c  Load the friction coefficients corresponding to side chains
 356       m=nct-nnt
 357       ind=0
 358       gamsc(21)=1.0d0
 359       do i=nnt,nct
 360         ind=ind+1
 361         ii = ind+m
 362         iti=itype(i)
 363         gamvec(ii)=gamsc(iti)
 364       enddo
 365       if (surfarea) call sdarea(gamvec)
 366 c      if (lprn) then
 367 c        write (iout,*) "Matrix A and vector gamma"
 368 c        do i=1,dimen1
 369 c          write (iout,'(i2,$)') i
 370 c          do j=1,dimen
 371 c            write (iout,'(f4.1,$)') A(i,j)
 372 c          enddo
 373 c          write (iout,'(f8.3)') gamvec(i)
 374 c        enddo
 375 c      endif
 376       if (lprn) then
 377         write (iout,*) "Vector gamvec"
 378         do i=1,dimen1
 379           write (iout,'(i5,f10.5)') i, gamvec(i)
 380         enddo
 381       endif
 382
 383 c The friction matrix
 384       do k=1,dimen
 385        do i=1,dimen
 386          dtdi=0.0d0
 387          do j=1,dimen1
 388            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
 389          enddo
 390          fricmat(k,i)=dtdi
 391        enddo
 392       enddo
 393
 394       if (lprn) then
 395         write (iout,'(//a)') "Matrix fricmat"
 396         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
 397       endif
 398       if (lang.eq.2 .or. lang.eq.3) then
 399 c Mass-scale the friction matrix if non-direct integration will be performed
 400       do i=1,dimen
 401         do j=1,dimen
 402           Ginvfric(i,j)=0.0d0
 403           do k=1,dimen
 404             do l=1,dimen
 405               Ginvfric(i,j)=Ginvfric(i,j)+
 406      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
 407             enddo
 408           enddo
 409         enddo
 410       enddo
 411 c Diagonalize the friction matrix
 412       ind=0
 413       do i=1,dimen
 414         do j=1,i
 415           ind=ind+1
 416           Ghalf(ind)=Ginvfric(i,j)
 417         enddo
 418       enddo
 419       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 420      &  ierr,iwork)
 421       if (lprn) then
 422         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 423      &    " mass-scaled friction matrix"
 424         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 425       endif
 426 c Precompute matrices for tinker stochastic integrator
 427 #ifndef LANG0
 428       do i=1,dimen
 429         do j=1,dimen
 430           mt1(i,j)=0.0d0
 431           mt2(i,j)=0.0d0
 432           do k=1,dimen
 433             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
 434             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
 435           enddo
 436           mt3(j,i)=mt1(i,j)
 437         enddo
 438       enddo
 439 #endif
 440       else if (lang.eq.4) then
 441 c Diagonalize the friction matrix
 442       ind=0
 443       do i=1,dimen
 444         do j=1,i
 445           ind=ind+1
 446           Ghalf(ind)=fricmat(i,j)
 447         enddo
 448       enddo
 449       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 450      &  ierr,iwork)
 451       if (lprn) then
 452         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 453      &    " friction matrix"
 454         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 455       endif
 456 c Determine the number of zero eigenvalues of the friction matrix
 457       nzero=max0(dimen-dimen1,0)
 458 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
 459 c        nzero=nzero+1
 460 c      enddo
 461       write (iout,*) "Number of zero eigenvalues:",nzero
 462       do i=1,dimen
 463         do j=1,dimen
 464           fricmat(i,j)=0.0d0
 465           do k=nzero+1,dimen
 466             fricmat(i,j)=fricmat(i,j)
 467      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k)
 468           enddo
 469         enddo
 470       enddo
 471       if (lprn) then
 472         write (iout,'(//a)') "Generalized inverse of fricmat"
 473         call matout(dimen,dimen,maxres6,maxres6,fricmat)
 474       endif
 475       endif
 476 #ifdef MPI
 477   10  continue
 478       if (nfgtasks.gt.1) then
 479         if (fg_rank.eq.0) then
 480 c The matching BROADCAST for fg processors is called in ERGASTULUM
 481           time00=MPI_Wtime()
 482           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
 483           time_Bcast=time_Bcast+MPI_Wtime()-time00
 484 c          print *,"Processor",myrank,
 485 c     &       " BROADCAST iorder in SETUP_FRICMAT"
 486         endif
 487 c       licznik=licznik+1
 488 c        write (iout,*) "setup_fricmat licznik",licznik
 489         time00=MPI_Wtime()
 490 c Scatter the friction matrix
 491         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
 492      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
 493      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
 494         time_scatter=time_scatter+MPI_Wtime()-time00
 495 #ifdef TIMING
 496         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
 497 #endif
 498         do i=1,dimen
 499           do j=1,2*my_ng_count
 500             fricmat(j,i)=fcopy(i,j)
 501           enddo
 502         enddo
 503 c        write (iout,*) "My chunk of fricmat"
 504 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
 505       endif
 506 #endif
 507       return
 508       end
 509 c-------------------------------------------------------------------------------
 510       subroutine sdarea(gamvec)
 511 c
 512 c Scale the friction coefficients according to solvent accessible surface areas
 513 c Code adapted from TINKER
 514 c AL 9/3/04
 515 c
 516       implicit real*8 (a-h,o-z)
 517       include 'DIMENSIONS'
 518       include 'COMMON.CONTROL'
 519       include 'COMMON.VAR'
 520       include 'COMMON.MD'
 521 #ifndef LANG0
 522       include 'COMMON.LANGEVIN'
 523 #else
 524       include 'COMMON.LANGEVIN.lang0'
 525 #endif
 526       include 'COMMON.CHAIN'
 527       include 'COMMON.DERIV'
 528       include 'COMMON.GEO'
 529       include 'COMMON.LOCAL'
 530       include 'COMMON.INTERACT'
 531       include 'COMMON.IOUNITS'
 532       include 'COMMON.NAMES'
 533       double precision radius(maxres2),gamvec(maxres6)
 534       parameter (twosix=1.122462048309372981d0)
 535       logical lprn /.true./
 536 c
 537 c     determine new friction coefficients every few SD steps
 538 c
 539 c     set the atomic radii to estimates of sigma values
 540 c
 541 c      print *,"Entered sdarea"
 542       probe = 0.0d0
 543
 544       do i=1,2*nres
 545         radius(i)=0.0d0
 546       enddo
 547 c  Load peptide group radii
 548       do i=nnt,nct-1
 549         radius(i)=pstok
 550       enddo
 551 c  Load side chain radii
 552       do i=nnt,nct
 553         iti=itype(i)
 554         radius(i+nres)=restok(iti)
 555       enddo
 556 c      do i=1,2*nres
 557 c        write (iout,*) "i",i," radius",radius(i)
 558 c      enddo
 559       do i = 1, 2*nres
 560          radius(i) = radius(i) / twosix
 561          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
 562       end do
 563 c
 564 c     scale atomic friction coefficients by accessible area
 565 c
 566       if (lprn) write (iout,*)
 567      &  "Original gammas, surface areas, scaling factors, new gammas, ",
 568      &  "std's of stochastic forces"
 569       ind=0
 570       do i=nnt,nct-1
 571         if (radius(i).gt.0.0d0) then
 572           call surfatom (i,area,radius)
 573           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
 574           if (lprn) write (iout,'(i5,3f10.5,$)')
 575      &      i,gamvec(ind+1),area,ratio
 576           do j=1,3
 577             ind=ind+1
 578             gamvec(ind) = ratio * gamvec(ind)
 579           enddo
 580           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
 581           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
 582         endif
 583       enddo
 584       do i=nnt,nct
 585         if (radius(i+nres).gt.0.0d0) then
 586           call surfatom (i+nres,area,radius)
 587           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
 588           if (lprn) write (iout,'(i5,3f10.5,$)')
 589      &      i,gamvec(ind+1),area,ratio
 590           do j=1,3
 591             ind=ind+1
 592             gamvec(ind) = ratio * gamvec(ind)
 593           enddo
 594           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
 595           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
 596         endif
 597       enddo
 598
 599       return
 600       end