source/unres/src_MD/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) 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) 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.21) 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) 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) 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       do i=nnt,nct
 369         ind=ind+1
 370         ii = ind+m
 371         iti=itype(i)
 372         gamvec(ii)=gamsc(iti)
 373       enddo
 374       if (surfarea) call sdarea(gamvec)
 375 c      if (lprn) then
 376 c        write (iout,*) "Matrix A and vector gamma"
 377 c        do i=1,dimen1
 378 c          write (iout,'(i2,$)') i
 379 c          do j=1,dimen
 380 c            write (iout,'(f4.1,$)') A(i,j)
 381 c          enddo
 382 c          write (iout,'(f8.3)') gamvec(i)
 383 c        enddo
 384 c      endif
 385       if (lprn) then
 386         write (iout,*) "Vector gamvec"
 387         do i=1,dimen1
 388           write (iout,'(i5,f10.5)') i, gamvec(i)
 389         enddo
 390       endif
 391
 392 c The friction matrix
 393       do k=1,dimen
 394        do i=1,dimen
 395          dtdi=0.0d0
 396          do j=1,dimen1
 397            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
 398          enddo
 399          fricmat(k,i)=dtdi
 400        enddo
 401       enddo
 402
 403       if (lprn) then
 404         write (iout,'(//a)') "Matrix fricmat"
 405         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
 406       endif
 407       if (lang.eq.2 .or. lang.eq.3) then
 408 c Mass-scale the friction matrix if non-direct integration will be performed
 409       do i=1,dimen
 410         do j=1,dimen
 411           Ginvfric(i,j)=0.0d0
 412           do k=1,dimen
 413             do l=1,dimen
 414               Ginvfric(i,j)=Ginvfric(i,j)+
 415      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
 416             enddo
 417           enddo
 418         enddo
 419       enddo
 420 c Diagonalize the friction matrix
 421       ind=0
 422       do i=1,dimen
 423         do j=1,i
 424           ind=ind+1
 425           Ghalf(ind)=Ginvfric(i,j)
 426         enddo
 427       enddo
 428       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 429      &  ierr,iwork)
 430       if (lprn) then
 431         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 432      &    " mass-scaled friction matrix"
 433         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 434       endif
 435 c Precompute matrices for tinker stochastic integrator
 436 #ifndef LANG0
 437       do i=1,dimen
 438         do j=1,dimen
 439           mt1(i,j)=0.0d0
 440           mt2(i,j)=0.0d0
 441           do k=1,dimen
 442             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
 443             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
 444           enddo
 445           mt3(j,i)=mt1(i,j)
 446         enddo
 447       enddo
 448 #endif
 449       else if (lang.eq.4) then
 450 c Diagonalize the friction matrix
 451       ind=0
 452       do i=1,dimen
 453         do j=1,i
 454           ind=ind+1
 455           Ghalf(ind)=fricmat(i,j)
 456         enddo
 457       enddo
 458       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 459      &  ierr,iwork)
 460       if (lprn) then
 461         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 462      &    " friction matrix"
 463         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 464       endif
 465 c Determine the number of zero eigenvalues of the friction matrix
 466       nzero=max0(dimen-dimen1,0)
 467 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
 468 c        nzero=nzero+1
 469 c      enddo
 470       write (iout,*) "Number of zero eigenvalues:",nzero
 471       do i=1,dimen
 472         do j=1,dimen
 473           fricmat(i,j)=0.0d0
 474           do k=nzero+1,dimen
 475             fricmat(i,j)=fricmat(i,j)
 476      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k)
 477           enddo
 478         enddo
 479       enddo
 480       if (lprn) then
 481         write (iout,'(//a)') "Generalized inverse of fricmat"
 482         call matout(dimen,dimen,maxres6,maxres6,fricmat)
 483       endif
 484       endif
 485 #ifdef MPI
 486   10  continue
 487       if (nfgtasks.gt.1) then
 488         if (fg_rank.eq.0) then
 489 c The matching BROADCAST for fg processors is called in ERGASTULUM
 490 #ifdef MPI
 491           time00=MPI_Wtime()
 492 #else
 493           time00=tcpu()
 494 #endif
 495           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
 496 #ifdef MPI
 497           time_Bcast=time_Bcast+MPI_Wtime()-time00
 498 #else
 499           time_Bcast=time_Bcast+tcpu()-time00
 500 #endif
 501 c          print *,"Processor",myrank,
 502 c     &       " BROADCAST iorder in SETUP_FRICMAT"
 503         endif
 504 c       licznik=licznik+1
 505 c        write (iout,*) "setup_fricmat licznik",licznik
 506 #ifdef MPI
 507         time00=MPI_Wtime()
 508 #else
 509         time00=tcpu()
 510 #endif
 511 c Scatter the friction matrix
 512         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
 513      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
 514      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
 515 #ifdef TIMING
 516 #ifdef MPI
 517         time_scatter=time_scatter+MPI_Wtime()-time00
 518         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
 519 #else
 520         time_scatter=time_scatter+tcpu()-time00
 521         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
 522 #endif
 523 #endif
 524         do i=1,dimen
 525           do j=1,2*my_ng_count
 526             fricmat(j,i)=fcopy(i,j)
 527           enddo
 528         enddo
 529 c        write (iout,*) "My chunk of fricmat"
 530 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
 531       endif
 532 #endif
 533       return
 534       end
 535 c-------------------------------------------------------------------------------
 536       subroutine sdarea(gamvec)
 537 c
 538 c Scale the friction coefficients according to solvent accessible surface areas
 539 c Code adapted from TINKER
 540 c AL 9/3/04
 541 c
 542       implicit real*8 (a-h,o-z)
 543       include 'DIMENSIONS'
 544       include 'COMMON.CONTROL'
 545       include 'COMMON.VAR'
 546       include 'COMMON.MD'
 547 #ifndef LANG0
 548       include 'COMMON.LANGEVIN'
 549 #else
 550       include 'COMMON.LANGEVIN.lang0'
 551 #endif
 552       include 'COMMON.CHAIN'
 553       include 'COMMON.DERIV'
 554       include 'COMMON.GEO'
 555       include 'COMMON.LOCAL'
 556       include 'COMMON.INTERACT'
 557       include 'COMMON.IOUNITS'
 558       include 'COMMON.NAMES'
 559       double precision radius(maxres2),gamvec(maxres6)
 560       parameter (twosix=1.122462048309372981d0)
 561       logical lprn /.false./
 562 c
 563 c     determine new friction coefficients every few SD steps
 564 c
 565 c     set the atomic radii to estimates of sigma values
 566 c
 567 c      print *,"Entered sdarea"
 568       probe = 0.0d0
 569
 570       do i=1,2*nres
 571         radius(i)=0.0d0
 572       enddo
 573 c  Load peptide group radii
 574       do i=nnt,nct-1
 575         radius(i)=pstok
 576       enddo
 577 c  Load side chain radii
 578       do i=nnt,nct
 579         iti=itype(i)
 580         radius(i+nres)=restok(iti)
 581       enddo
 582 c      do i=1,2*nres
 583 c        write (iout,*) "i",i," radius",radius(i)
 584 c      enddo
 585       do i = 1, 2*nres
 586          radius(i) = radius(i) / twosix
 587          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
 588       end do
 589 c
 590 c     scale atomic friction coefficients by accessible area
 591 c
 592       if (lprn) write (iout,*)
 593      &  "Original gammas, surface areas, scaling factors, new gammas, ",
 594      &  "std's of stochastic forces"
 595       ind=0
 596       do i=nnt,nct-1
 597         if (radius(i).gt.0.0d0) then
 598           call surfatom (i,area,radius)
 599           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
 600           if (lprn) write (iout,'(i5,3f10.5,$)')
 601      &      i,gamvec(ind+1),area,ratio
 602           do j=1,3
 603             ind=ind+1
 604             gamvec(ind) = ratio * gamvec(ind)
 605           enddo
 606           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
 607           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
 608         endif
 609       enddo
 610       do i=nnt,nct
 611         if (radius(i+nres).gt.0.0d0) then
 612           call surfatom (i+nres,area,radius)
 613           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
 614           if (lprn) write (iout,'(i5,3f10.5,$)')
 615      &      i,gamvec(ind+1),area,ratio
 616           do j=1,3
 617             ind=ind+1
 618             gamvec(ind) = ratio * gamvec(ind)
 619           enddo
 620           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
 621           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
 622         endif
 623       enddo
 624
 625       return
 626       end