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