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.ntyp1)) 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.ntyp1)) 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 #ifdef MPI
 189       time00=MPI_Wtime()
 190 #else
 191       time00=tcpu()
 192 #endif
 193 c Compute the stochastic forces acting on bodies. Store in force.
 194       do i=nnt,nct-1
 195         sig=stdforcp(i)
 196         lowb=-5*sig
 197         highb=5*sig
 198         do j=1,3
 199           force(j,i)=anorm_distr(x,sig,lowb,highb)
 200         enddo
 201       enddo
 202       do i=nnt,nct
 203         sig2=stdforcsc(i)
 204         lowb2=-5*sig2
 205         highb2=5*sig2
 206         do j=1,3
 207           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
 208         enddo
 209       enddo
 210 #ifdef MPI
 211       time_fsample=time_fsample+MPI_Wtime()-time00
 212 #else
 213       time_fsample=time_fsample+tcpu()-time00
 214 #endif
 215 c Compute the stochastic forces acting on virtual-bond vectors.
 216       do j=1,3
 217         ff(j)=0.0d0
 218       enddo
 219       do i=nct-1,nnt,-1
 220         do j=1,3
 221           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
 222         enddo
 223         do j=1,3
 224           ff(j)=ff(j)+force(j,i)
 225         enddo
 226         if (itype(i+1).ne.ntyp1) then
 227           do j=1,3
 228             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
 229             ff(j)=ff(j)+force(j,i+nres+1)
 230           enddo
 231         endif
 232       enddo
 233       do j=1,3
 234         stochforc(j,0)=ff(j)+force(j,nnt+nres)
 235       enddo
 236       do i=nnt,nct
 237         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 238           do j=1,3
 239             stochforc(j,i+nres)=force(j,i+nres)
 240           enddo
 241         endif
 242       enddo
 243
 244       do j=1,3
 245         stochforcvec(j)=stochforc(j,0)
 246       enddo
 247       ind=3
 248       do i=nnt,nct-1
 249         do j=1,3
 250           stochforcvec(ind+j)=stochforc(j,i)
 251         enddo
 252         ind=ind+3
 253       enddo
 254       do i=nnt,nct
 255         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 256           do j=1,3
 257             stochforcvec(ind+j)=stochforc(j,i+nres)
 258           enddo
 259           ind=ind+3
 260         endif
 261       enddo
 262       if (lprn) then
 263         write (iout,*) "stochforcvec"
 264         do i=1,3*dimen
 265           write(iout,'(i5,e15.5)') i,stochforcvec(i)
 266         enddo
 267         write(iout,*) "Stochastic forces backbone"
 268         do i=0,nct-1
 269           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
 270         enddo
 271         write(iout,*) "Stochastic forces side chain"
 272         do i=nnt,nct
 273           write(iout,'(i5,3e15.5)')
 274      &      i,(stochforc(j,i+nres),j=1,3)
 275         enddo
 276       endif
 277
 278       if (lprn) then
 279
 280       ind=0
 281       do i=nnt,nct-1
 282         write (iout,*) i,ind
 283         do j=1,3
 284           forcvec(ind+j)=force(j,i)
 285         enddo
 286         ind=ind+3
 287       enddo
 288       do i=nnt,nct
 289         write (iout,*) i,ind
 290         do j=1,3
 291           forcvec(j+ind)=force(j,i+nres)
 292         enddo
 293         ind=ind+3
 294       enddo
 295
 296       write (iout,*) "forcvec"
 297       ind=0
 298       do i=nnt,nct-1
 299         do j=1,3
 300           write (iout,'(2i3,2f10.5)') i,j,force(j,i),
 301      &      forcvec(ind+j)
 302         enddo
 303         ind=ind+3
 304       enddo
 305       do i=nnt,nct
 306         do j=1,3
 307           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
 308      &     forcvec(ind+j)
 309         enddo
 310         ind=ind+3
 311       enddo
 312
 313       endif
 314
 315       return
 316       end
 317 c------------------------------------------------------------------
 318       subroutine setup_fricmat
 319       implicit real*8 (a-h,o-z)
 320 #ifdef MPI
 321       include 'mpif.h'
 322 #endif
 323       include 'DIMENSIONS'
 324       include 'COMMON.VAR'
 325       include 'COMMON.CHAIN'
 326       include 'COMMON.DERIV'
 327       include 'COMMON.GEO'
 328       include 'COMMON.LOCAL'
 329       include 'COMMON.INTERACT'
 330       include 'COMMON.MD'
 331       include 'COMMON.SETUP'
 332       include 'COMMON.TIME1'
 333 c      integer licznik /0/
 334 c      save licznik
 335 #ifndef LANG0
 336       include 'COMMON.LANGEVIN'
 337 #else
 338       include 'COMMON.LANGEVIN.lang0'
 339 #endif
 340       include 'COMMON.IOUNITS'
 341       integer IERROR
 342       integer i,j,ind,ind1,m
 343       logical lprn /.false./
 344       double precision dtdi,gamvec(MAXRES2),
 345      &  ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2)
 346       common /syfek/ gamvec
 347       double precision work(8*maxres2)
 348       integer iwork(maxres2)
 349       common /przechowalnia/ ginvfric,Ghalf,fcopy
 350 #ifdef MPI
 351       if (fg_rank.ne.king) goto 10
 352 #endif
 353 c  Zeroing out fricmat
 354       do i=1,dimen
 355         do j=1,dimen
 356           fricmat(i,j)=0.0d0
 357         enddo
 358       enddo
 359 c  Load the friction coefficients corresponding to peptide groups
 360       ind1=0
 361       do i=nnt,nct-1
 362         ind1=ind1+1
 363         gamvec(ind1)=gamp
 364       enddo
 365 c  Load the friction coefficients corresponding to side chains
 366       m=nct-nnt
 367       ind=0
 368 C      gamsc(ntyp1)=1.0d0
 369       do i=nnt,nct
 370         ind=ind+1
 371         ii = ind+m
 372         iti=itype(i)
 373         gamvec(ii)=gamsc(iabs(iti))
 374       enddo
 375       if (surfarea) call sdarea(gamvec)
 376 c      if (lprn) then
 377 c        write (iout,*) "Matrix A and vector gamma"
 378 c        do i=1,dimen1
 379 c          write (iout,'(i2,$)') i
 380 c          do j=1,dimen
 381 c            write (iout,'(f4.1,$)') A(i,j)
 382 c          enddo
 383 c          write (iout,'(f8.3)') gamvec(i)
 384 c        enddo
 385 c      endif
 386       if (lprn) then
 387         write (iout,*) "Vector gamvec"
 388         do i=1,dimen1
 389           write (iout,'(i5,f10.5)') i, gamvec(i)
 390         enddo
 391       endif
 392
 393 c The friction matrix
 394       do k=1,dimen
 395        do i=1,dimen
 396          dtdi=0.0d0
 397          do j=1,dimen1
 398            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
 399          enddo
 400          fricmat(k,i)=dtdi
 401        enddo
 402       enddo
 403
 404       if (lprn) then
 405         write (iout,'(//a)') "Matrix fricmat"
 406         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
 407       endif
 408       if (lang.eq.2 .or. lang.eq.3) then
 409 c Mass-scale the friction matrix if non-direct integration will be performed
 410       do i=1,dimen
 411         do j=1,dimen
 412           Ginvfric(i,j)=0.0d0
 413           do k=1,dimen
 414             do l=1,dimen
 415               Ginvfric(i,j)=Ginvfric(i,j)+
 416      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
 417             enddo
 418           enddo
 419         enddo
 420       enddo
 421 c Diagonalize the friction matrix
 422       ind=0
 423       do i=1,dimen
 424         do j=1,i
 425           ind=ind+1
 426           Ghalf(ind)=Ginvfric(i,j)
 427         enddo
 428       enddo
 429       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 430      &  ierr,iwork)
 431       if (lprn) then
 432         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 433      &    " mass-scaled friction matrix"
 434         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 435       endif
 436 c Precompute matrices for tinker stochastic integrator
 437 #ifndef LANG0
 438       do i=1,dimen
 439         do j=1,dimen
 440           mt1(i,j)=0.0d0
 441           mt2(i,j)=0.0d0
 442           do k=1,dimen
 443             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
 444             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
 445           enddo
 446           mt3(j,i)=mt1(i,j)
 447         enddo
 448       enddo
 449 #endif
 450       else if (lang.eq.4) then
 451 c Diagonalize the friction matrix
 452       ind=0
 453       do i=1,dimen
 454         do j=1,i
 455           ind=ind+1
 456           Ghalf(ind)=fricmat(i,j)
 457         enddo
 458       enddo
 459       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 460      &  ierr,iwork)
 461       if (lprn) then
 462         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 463      &    " friction matrix"
 464         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 465       endif
 466 c Determine the number of zero eigenvalues of the friction matrix
 467       nzero=max0(dimen-dimen1,0)
 468 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
 469 c        nzero=nzero+1
 470 c      enddo
 471       write (iout,*) "Number of zero eigenvalues:",nzero
 472       do i=1,dimen
 473         do j=1,dimen
 474           fricmat(i,j)=0.0d0
 475           do k=nzero+1,dimen
 476             fricmat(i,j)=fricmat(i,j)
 477      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k)
 478           enddo
 479         enddo
 480       enddo
 481       if (lprn) then
 482         write (iout,'(//a)') "Generalized inverse of fricmat"
 483         call matout(dimen,dimen,maxres6,maxres6,fricmat)
 484       endif
 485       endif
 486 #ifdef MPI
 487   10  continue
 488       if (nfgtasks.gt.1) then
 489         if (fg_rank.eq.0) then
 490 c The matching BROADCAST for fg processors is called in ERGASTULUM
 491 #ifdef MPI
 492           time00=MPI_Wtime()
 493 #else
 494           time00=tcpu()
 495 #endif
 496           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
 497 #ifdef MPI
 498           time_Bcast=time_Bcast+MPI_Wtime()-time00
 499 #else
 500           time_Bcast=time_Bcast+tcpu()-time00
 501 #endif
 502 c          print *,"Processor",myrank,
 503 c     &       " BROADCAST iorder in SETUP_FRICMAT"
 504         endif
 505 c       licznik=licznik+1
 506 c        write (iout,*) "setup_fricmat licznik",licznik
 507 #ifdef MPI
 508         time00=MPI_Wtime()
 509 #else
 510         time00=tcpu()
 511 #endif
 512 c Scatter the friction matrix
 513         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
 514      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
 515      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
 516 #ifdef TIMING
 517 #ifdef MPI
 518         time_scatter=time_scatter+MPI_Wtime()-time00
 519         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
 520 #else
 521         time_scatter=time_scatter+tcpu()-time00
 522         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
 523 #endif
 524 #endif
 525         do i=1,dimen
 526           do j=1,2*my_ng_count
 527             fricmat(j,i)=fcopy(i,j)
 528           enddo
 529         enddo
 530 c        write (iout,*) "My chunk of fricmat"
 531 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
 532       endif
 533 #endif
 534       return
 535       end
 536 c-------------------------------------------------------------------------------
 537       subroutine sdarea(gamvec)
 538 c
 539 c Scale the friction coefficients according to solvent accessible surface areas
 540 c Code adapted from TINKER
 541 c AL 9/3/04
 542 c
 543       implicit real*8 (a-h,o-z)
 544       include 'DIMENSIONS'
 545       include 'COMMON.CONTROL'
 546       include 'COMMON.VAR'
 547       include 'COMMON.MD'
 548 #ifndef LANG0
 549       include 'COMMON.LANGEVIN'
 550 #else
 551       include 'COMMON.LANGEVIN.lang0'
 552 #endif
 553       include 'COMMON.CHAIN'
 554       include 'COMMON.DERIV'
 555       include 'COMMON.GEO'
 556       include 'COMMON.LOCAL'
 557       include 'COMMON.INTERACT'
 558       include 'COMMON.IOUNITS'
 559       include 'COMMON.NAMES'
 560       double precision radius(maxres2),gamvec(maxres2)
 561       parameter (twosix=1.122462048309372981d0)
 562       logical lprn /.false./
 563 c
 564 c     determine new friction coefficients every few SD steps
 565 c
 566 c     set the atomic radii to estimates of sigma values
 567 c
 568 c      print *,"Entered sdarea"
 569       probe = 0.0d0
 570
 571       do i=1,2*nres
 572         radius(i)=0.0d0
 573       enddo
 574 c  Load peptide group radii
 575       do i=nnt,nct-1
 576         radius(i)=pstok
 577       enddo
 578 c  Load side chain radii
 579       do i=nnt,nct
 580         iti=itype(i)
 581         radius(i+nres)=restok(iti)
 582       enddo
 583 c      do i=1,2*nres
 584 c        write (iout,*) "i",i," radius",radius(i)
 585 c      enddo
 586       do i = 1, 2*nres
 587          radius(i) = radius(i) / twosix
 588          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
 589       end do
 590 c
 591 c     scale atomic friction coefficients by accessible area
 592 c
 593       if (lprn) write (iout,*)
 594      &  "Original gammas, surface areas, scaling factors, new gammas, ",
 595      &  "std's of stochastic forces"
 596       ind=0
 597       do i=nnt,nct-1
 598         if (radius(i).gt.0.0d0) then
 599           call surfatom (i,area,radius)
 600           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
 601           if (lprn) write (iout,'(i5,3f10.5,$)')
 602      &      i,gamvec(ind+1),area,ratio
 603           do j=1,3
 604             ind=ind+1
 605             gamvec(ind) = ratio * gamvec(ind)
 606           enddo
 607           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
 608           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
 609         endif
 610       enddo
 611       do i=nnt,nct
 612         if (radius(i+nres).gt.0.0d0) then
 613           call surfatom (i+nres,area,radius)
 614           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
 615           if (lprn) write (iout,'(i5,3f10.5,$)')
 616      &      i,gamvec(ind+1),area,ratio
 617           do j=1,3
 618             ind=ind+1
 619             gamvec(ind) = ratio * gamvec(ind)
 620           enddo
 621           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
 622           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
 623         endif
 624       enddo
 625
 626       return
 627       end