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