source/unres/src-HCD-5D/stochfric.F

   1       subroutine friction_force
   2       implicit none
   3 #ifdef MPI
   4       include 'mpif.h'
   5 #endif
   6       include 'DIMENSIONS'
   7       include 'COMMON.VAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.DERIV'
  10       include 'COMMON.GEO'
  11       include 'COMMON.LOCAL'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.MD'
  14 #ifdef FIVEDIAG
  15        include 'COMMON.LAGRANGE.5diag'
  16 #else
  17        include 'COMMON.LAGRANGE'
  18 #endif
  19 #ifndef LANG0
  20       include 'COMMON.LANGEVIN'
  21 #else
  22 #ifdef FIVEDIAG
  23       include 'COMMON.LANGEVIN.lang0.5diag'
  24 #else
  25       include 'COMMON.LANGEVIN.lang0'
  26 #endif
  27 #endif
  28       include 'COMMON.IOUNITS'
  29 #ifdef FIVEDIAG
  30       integer iposc,ichain,n,innt,inct
  31       double precision v_work(3,maxres2),vvec(maxres2),rs(maxres2)
  32 #else
  33       double precision gamvec(MAXRES6)
  34       common /syfek/ gamvec
  35       double precision vv(3),vvtot(3,maxres),v_work(MAXRES6),
  36      & ginvfric(maxres2,maxres2)
  37       common /przechowalnia/ ginvfric
  38 #endif
  39       integer i,j,k,ind
  40
  41       logical lprn /.false./, checkmode /.false./
  42 #ifdef FIVEDIAG
  43 #ifdef TIMING
  44       include 'COMMON.TIME1'
  45       double precision time01
  46 #endif
  47 c Here accelerations due to friction forces are computed right after forces.
  48       d_t_work(:6*nres)=0.0d0
  49       do j=1,3
  50         v_work(j,1)=d_t(j,0)
  51         v_work(j,nnt)=d_t(j,0)
  52       enddo
  53       do i=nnt+1,nct
  54         do j=1,3
  55           v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
  56         enddo
  57       enddo
  58       do i=nnt,nct
  59         if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
  60           do j=1,3
  61             v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
  62           enddo
  63         endif
  64       enddo
  65 #ifdef DEBUG
  66       write (iout,*) "v_work"
  67       do i=1,2*nres
  68         write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
  69       enddo
  70 #endif
  71       do j=1,3
  72         ind=0
  73         do ichain=1,nchain
  74           n=dimen_chain(ichain)
  75           iposc=iposd_chain(ichain)
  76 c          write (iout,*) "friction_force j",j," ichain",ichain,
  77 c     &       " n",n," iposc",iposc,iposc+n-1
  78           innt=chain_border(1,ichain)
  79           inct=chain_border(2,ichain)
  80 c diagnostics
  81 c          innt=chain_border(1,1)
  82 c          inct=chain_border(2,1)
  83           do i=innt,inct
  84             vvec(ind+1)=v_work(j,i)
  85             ind=ind+1
  86             if (iabs(itype(i)).ne.10) then
  87               vvec(ind+1)=v_work(j,i+nres)
  88               ind=ind+1
  89             endif
  90           enddo
  91 #ifdef DEBUG
  92           write (iout,*) "vvec ind",ind," n",n
  93           write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
  94 #endif
  95 c          write (iout,*) "chain",i," ind",ind," n",n
  96 #ifdef TIMING
  97           time01=MPI_Wtime()
  98 #endif
  99           call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
 100      &     DU2fric(iposc),vvec(iposc),rs)
 101 #ifdef TIMING
 102           time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 103 #endif
 104 #ifdef DEBUG
 105           write (iout,*) "rs"
 106           write (iout,'(f10.5)') (rs(i),i=1,n)
 107 #endif
 108           do i=iposc,iposc+n-1
 109 c            write (iout,*) "ichain",ichain," i",i," j",j,
 110 c     &       "index",3*(i-1)+j,"rs",rs(i-iposc+1)
 111             fric_work(3*(i-1)+j)=-rs(i-iposc+1)
 112           enddo
 113         enddo
 114       enddo
 115 #ifdef DEBUG
 116       write (iout,*) "Vector fric_work dimen3",dimen3
 117       write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
 118 #endif
 119 #else
 120       do i=0,2*nres
 121         do j=1,3
 122           friction(j,i)=0.0d0
 123         enddo
 124       enddo
 125
 126       do j=1,3
 127         d_t_work(j)=d_t(j,0)
 128       enddo
 129       ind=3
 130       do i=nnt,nct-1
 131         do j=1,3
 132           d_t_work(ind+j)=d_t(j,i)
 133         enddo
 134         ind=ind+3
 135       enddo
 136       do i=nnt,nct
 137         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 138           do j=1,3
 139             d_t_work(ind+j)=d_t(j,i+nres)
 140           enddo
 141           ind=ind+3
 142         endif
 143       enddo
 144
 145       call fricmat_mult(d_t_work,fric_work)
 146
 147       if (.not.checkmode) return
 148
 149       if (lprn) then
 150         write (iout,*) "d_t_work and fric_work"
 151         do i=1,3*dimen
 152           write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
 153         enddo
 154       endif
 155       do j=1,3
 156         friction(j,0)=fric_work(j)
 157       enddo
 158       ind=3
 159       do i=nnt,nct-1
 160         do j=1,3
 161           friction(j,i)=fric_work(ind+j)
 162         enddo
 163         ind=ind+3
 164       enddo
 165       do i=nnt,nct
 166         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 167           do j=1,3
 168             friction(j,i+nres)=fric_work(ind+j)
 169           enddo
 170           ind=ind+3
 171         endif
 172       enddo
 173       if (lprn) then
 174         write(iout,*) "Friction backbone"
 175         do i=0,nct-1
 176           write(iout,'(i5,3e15.5,5x,3e15.5)')
 177      &     i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
 178         enddo
 179         write(iout,*) "Friction side chain"
 180         do i=nnt,nct
 181           write(iout,'(i5,3e15.5,5x,3e15.5)')
 182      &     i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
 183         enddo
 184       endif
 185       if (lprn) then
 186         do j=1,3
 187           vv(j)=d_t(j,0)
 188         enddo
 189         do i=nnt,nct
 190           do j=1,3
 191             vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
 192             vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
 193             vv(j)=vv(j)+d_t(j,i)
 194           enddo
 195         enddo
 196         write (iout,*) "vvtot backbone and sidechain"
 197         do i=nnt,nct
 198           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),
 199      &     (vvtot(j,i+nres),j=1,3)
 200         enddo
 201         ind=0
 202         do i=nnt,nct-1
 203           do j=1,3
 204             v_work(ind+j)=vvtot(j,i)
 205           enddo
 206           ind=ind+3
 207         enddo
 208         do i=nnt,nct
 209           do j=1,3
 210             v_work(ind+j)=vvtot(j,i+nres)
 211           enddo
 212           ind=ind+3
 213         enddo
 214         write (iout,*) "v_work gamvec and site-based friction forces"
 215         do i=1,dimen1
 216           write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),
 217      &      gamvec(i)*v_work(i)
 218         enddo
 219 c        do i=1,dimen
 220 c          fric_work1(i)=0.0d0
 221 c          do j=1,dimen1
 222 c            fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
 223 c          enddo
 224 c        enddo
 225 c        write (iout,*) "fric_work and fric_work1"
 226 c        do i=1,dimen
 227 c          write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
 228 c        enddo
 229         do i=1,dimen
 230           do j=1,dimen
 231             ginvfric(i,j)=0.0d0
 232             do k=1,dimen
 233               ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
 234             enddo
 235           enddo
 236         enddo
 237         write (iout,*) "ginvfric"
 238         do i=1,dimen
 239           write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
 240         enddo
 241         write (iout,*) "symmetry check"
 242         do i=1,dimen
 243           do j=1,i-1
 244             write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
 245           enddo
 246         enddo
 247       endif
 248 #endif
 249       return
 250       end
 251 c-----------------------------------------------------
 252       subroutine stochastic_force(stochforcvec)
 253       implicit none
 254       include 'DIMENSIONS'
 255 #ifdef MPI
 256       include 'mpif.h'
 257       double precision time00
 258 #endif
 259       include 'COMMON.VAR'
 260       include 'COMMON.CHAIN'
 261       include 'COMMON.DERIV'
 262       include 'COMMON.GEO'
 263       include 'COMMON.LOCAL'
 264       include 'COMMON.INTERACT'
 265       include 'COMMON.MD'
 266 #ifdef FIVEDIAG
 267        include 'COMMON.LAGRANGE.5diag'
 268 #else
 269        include 'COMMON.LAGRANGE'
 270 #endif
 271       include 'COMMON.TIME1'
 272 #ifndef LANG0
 273       include 'COMMON.LANGEVIN'
 274 #else
 275 #ifdef FIVEDIAG
 276       include 'COMMON.LANGEVIN.lang0.5diag'
 277 #else
 278       include 'COMMON.LANGEVIN.lang0'
 279 #endif
 280 #endif
 281       include 'COMMON.IOUNITS'
 282
 283       double precision x,sig,lowb,highb,
 284      & ff(3),force(3,0:MAXRES2),zeta2,lowb2,
 285      & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6)
 286       logical lprn /.false./
 287       integer i,j,ind,ii,iti
 288       double precision anorm_distr
 289 #ifdef FIVEDIAG
 290       integer ichain,innt,inct,iposc
 291 #endif
 292
 293       do i=0,2*nres
 294         do j=1,3
 295           stochforc(j,i)=0.0d0
 296         enddo
 297       enddo
 298       x=0.0d0
 299
 300 #ifdef MPI
 301       time00=MPI_Wtime()
 302 #else
 303       time00=tcpu()
 304 #endif
 305 c Compute the stochastic forces acting on bodies. Store in force.
 306       do i=nnt,nct-1
 307 #ifdef FIVEDIAG
 308         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
 309 #endif
 310         sig=stdforcp(i)
 311         lowb=-5*sig
 312         highb=5*sig
 313         do j=1,3
 314           force(j,i)=anorm_distr(x,sig,lowb,highb)
 315         enddo
 316       enddo
 317       do i=nnt,nct
 318         sig2=stdforcsc(i)
 319         lowb2=-5*sig2
 320         highb2=5*sig2
 321         do j=1,3
 322           force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
 323         enddo
 324       enddo
 325 #ifdef DEBUG
 326       write (iout,*) "Stochastic forces on sites"
 327       do i=1,nres
 328         write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3),
 329      &     (force(j,i+nres),j=1,3)
 330       enddo
 331 #endif
 332 #ifdef MPI
 333       time_fsample=time_fsample+MPI_Wtime()-time00
 334 #else
 335       time_fsample=time_fsample+tcpu()-time00
 336 #endif
 337 #ifdef FIVEDIAG
 338       ind=0
 339       do ichain=1,nchain
 340         innt=chain_border(1,ichain)
 341         inct=chain_border(2,ichain)
 342         iposc=iposd_chain(ichain)
 343 c for debugging only
 344 c        innt=chain_border(1,1)
 345 c        inct=chain_border(2,1)
 346 c        iposc=iposd_chain(1)
 347 c        write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
 348 c     &    " inct",inct," iposc",iposc
 349         do j=1,3
 350           stochforcvec(ind+j)=0.5d0*force(j,innt)
 351         enddo
 352         if (iabs(itype(innt)).eq.10) then
 353           do j=1,3
 354             stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
 355           enddo
 356           ind=ind+3
 357         else
 358           ind=ind+3
 359           do j=1,3
 360             stochforcvec(ind+j)=force(j,innt+nres)
 361           enddo
 362           ind=ind+3
 363         endif
 364         do i=innt+1,inct-1
 365           do j=1,3
 366             stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
 367           enddo
 368           if (iabs(itype(i)).eq.10) then
 369             do j=1,3
 370               stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
 371             enddo
 372             ind=ind+3
 373           else
 374             ind=ind+3
 375             do j=1,3
 376               stochforcvec(ind+j)=force(j,i+nres)
 377             enddo
 378             ind=ind+3
 379           endif
 380         enddo
 381         do j=1,3
 382           stochforcvec(ind+j)=0.5d0*force(j,inct-1)
 383         enddo
 384         if (iabs(itype(inct)).eq.10) then
 385           do j=1,3
 386             stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
 387           enddo
 388           ind=ind+3
 389         else
 390           ind=ind+3
 391           do j=1,3
 392             stochforcvec(ind+j)=force(j,inct+nres)
 393           enddo
 394           ind=ind+3
 395         endif
 396 c        write (iout,*) "chain",ichain," ind",ind
 397       enddo
 398 #ifdef DEBUG
 399       write (iout,*) "stochforcvec"
 400       write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
 401 #endif
 402 #else
 403 c Compute the stochastic forces acting on virtual-bond vectors.
 404       do j=1,3
 405         ff(j)=0.0d0
 406       enddo
 407       do i=nct-1,nnt,-1
 408         do j=1,3
 409           stochforc(j,i)=ff(j)+0.5d0*force(j,i)
 410         enddo
 411         do j=1,3
 412           ff(j)=ff(j)+force(j,i)
 413         enddo
 414         if (itype(i+1).ne.ntyp1) then
 415           do j=1,3
 416             stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
 417             ff(j)=ff(j)+force(j,i+nres+1)
 418           enddo
 419         endif
 420       enddo
 421       do j=1,3
 422         stochforc(j,0)=ff(j)+force(j,nnt+nres)
 423       enddo
 424       do i=nnt,nct
 425         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 426           do j=1,3
 427             stochforc(j,i+nres)=force(j,i+nres)
 428           enddo
 429         endif
 430       enddo
 431       do j=1,3
 432         stochforcvec(j)=stochforc(j,0)
 433       enddo
 434       ind=3
 435       do i=nnt,nct-1
 436         do j=1,3
 437           stochforcvec(ind+j)=stochforc(j,i)
 438         enddo
 439         ind=ind+3
 440       enddo
 441       do i=nnt,nct
 442         if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 443           do j=1,3
 444             stochforcvec(ind+j)=stochforc(j,i+nres)
 445           enddo
 446           ind=ind+3
 447         endif
 448       enddo
 449       if (lprn) then
 450         write (iout,*) "stochforcvec"
 451         do i=1,3*dimen
 452           write(iout,'(i5,e15.5)') i,stochforcvec(i)
 453         enddo
 454         write(iout,*) "Stochastic forces backbone"
 455         do i=0,nct-1
 456           write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
 457         enddo
 458         write(iout,*) "Stochastic forces side chain"
 459         do i=nnt,nct
 460           write(iout,'(i5,3e15.5)')
 461      &      i,(stochforc(j,i+nres),j=1,3)
 462         enddo
 463       endif
 464
 465       if (lprn) then
 466
 467       ind=0
 468       do i=nnt,nct-1
 469         write (iout,*) i,ind
 470         do j=1,3
 471           forcvec(ind+j)=force(j,i)
 472         enddo
 473         ind=ind+3
 474       enddo
 475       do i=nnt,nct
 476         write (iout,*) i,ind
 477         do j=1,3
 478           forcvec(j+ind)=force(j,i+nres)
 479         enddo
 480         ind=ind+3
 481       enddo
 482
 483       write (iout,*) "forcvec"
 484       ind=0
 485       do i=nnt,nct-1
 486         do j=1,3
 487           write (iout,'(2i3,2f10.5)') i,j,force(j,i),
 488      &      forcvec(ind+j)
 489         enddo
 490         ind=ind+3
 491       enddo
 492       do i=nnt,nct
 493         do j=1,3
 494           write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),
 495      &     forcvec(ind+j)
 496         enddo
 497         ind=ind+3
 498       enddo
 499       endif
 500 #endif
 501       return
 502       end
 503 c------------------------------------------------------------------
 504       subroutine setup_fricmat
 505       implicit none
 506 #ifdef MPI
 507       include 'mpif.h'
 508       integer ierr
 509       double precision time00
 510 #endif
 511       include 'DIMENSIONS'
 512 #ifndef FIVEDIAG
 513       integer mmaxres2
 514       parameter (mmaxres2=(maxres2*(maxres2+1)/2))
 515 #endif
 516       include 'COMMON.VAR'
 517       include 'COMMON.CHAIN'
 518       include 'COMMON.DERIV'
 519       include 'COMMON.GEO'
 520       include 'COMMON.LOCAL'
 521       include 'COMMON.INTERACT'
 522       include 'COMMON.MD'
 523 #ifdef FIVEDIAG
 524        include 'COMMON.LAGRANGE.5diag'
 525 #else
 526        include 'COMMON.LAGRANGE'
 527 #endif
 528       include 'COMMON.SETUP'
 529       include 'COMMON.TIME1'
 530 c      integer licznik /0/
 531 c      save licznik
 532 #ifndef LANG0
 533       include 'COMMON.LANGEVIN'
 534 #else
 535 #ifdef FIVEDIAG
 536       include 'COMMON.LANGEVIN.lang0.5diag'
 537 #else
 538       include 'COMMON.LANGEVIN.lang0'
 539 #endif
 540 #endif
 541       include 'COMMON.IOUNITS'
 542       integer IERROR
 543       integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
 544       integer ichain,nind
 545       logical lprn /.false./
 546       double precision dtdi,gamvec(MAXRES2)
 547       common /syfek/ gamvec
 548 #ifndef FIVEDIAG
 549       double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
 550      & fcopy(maxres2,maxres2)
 551       double precision work(8*maxres2)
 552       integer iwork(maxres2)
 553       common /przechowalnia/ ginvfric,Ghalf,fcopy
 554 #endif
 555
 556 #ifdef MPI
 557       if (fg_rank.ne.king) goto 10
 558 #endif
 559       ind1=0
 560       do i=nnt,nct-1
 561         ind1=ind1+1
 562         gamvec(ind1)=gamp
 563       enddo
 564 c  Load the friction coefficients corresponding to side chains
 565       m=nct-nnt
 566       if (lprn) write (iout,*) "m",m
 567       ind=0
 568 C      gamsc(ntyp1)=1.0d0
 569       do i=nnt,nct
 570         ind=ind+1
 571         ii = ind+m
 572         iti=itype(i)
 573         gamvec(ii)=gamsc(iabs(iti))
 574       enddo
 575       if (surfarea) call sdarea(gamvec)
 576       if (lprn) then
 577         write (iout,*) "Vector gamvec ii",ii
 578         do i=1,ii
 579           write (iout,'(i5,f10.5)') i, gamvec(i)
 580         enddo
 581       endif
 582 #ifdef FIVEDIAG
 583       DMfric(:2*nres)=0.0d0
 584       DU1fric(:2*nres)=0.0d0
 585       DU2fric(:2*nres)=0.0d0
 586       ind=1
 587       do ichain=1,nchain
 588         innt=chain_border(1,ichain)
 589         inct=chain_border(2,ichain)
 590 c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
 591 c DMfric part
 592         DMfric(ind)=gamvec(innt-nnt+1)/4
 593         if (iabs(itype(innt)).eq.10) then
 594           DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
 595           ind=ind+1
 596         else
 597           DMfric(ind+1)=gamvec(m+innt-nnt+1)
 598           ind=ind+2
 599         endif
 600 c        write (iout,*) "DMfric init ind",ind
 601 c DMfric
 602         do i=innt+1,inct-1
 603           DMfric(ind)=gamvec(i-nnt+1)/2
 604           if (iabs(itype(i)).eq.10) then
 605             DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
 606             ind=ind+1
 607           else
 608             DMfric(ind+1)=gamvec(m+i-nnt+1)
 609             ind=ind+2
 610           endif
 611         enddo
 612 c        write (iout,*) "DMfric endloop ind",ind
 613         if (inct.gt.innt) then
 614           DMfric(ind)=gamvec(inct-1-nnt+1)/4
 615           if (iabs(itype(inct)).eq.10) then
 616             DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
 617             ind=ind+1
 618           else
 619             DMfric(ind+1)=gamvec(inct+m-nnt+1)
 620             ind=ind+2
 621           endif
 622         endif
 623 c        write (iout,*) "DMfric end ind",ind
 624       enddo
 625 c DU1fric part
 626       do ichain=1,nchain
 627         ind=iposd_chain(ichain)
 628         innt=chain_border(1,ichain)
 629         inct=chain_border(2,ichain)
 630         do i=innt,inct
 631           if (iabs(itype(i)).ne.10) then
 632             ind=ind+2
 633           else
 634             DU1fric(ind)=gamvec(i-nnt+1)/4
 635             ind=ind+1
 636           endif
 637         enddo
 638       enddo
 639 c DU2fric part
 640       do ichain=1,nchain
 641         ind=iposd_chain(ichain)
 642         innt=chain_border(1,ichain)
 643         inct=chain_border(2,ichain)
 644         do i=innt,inct-1
 645           if (iabs(itype(i)).ne.10) then
 646             DU2fric(ind)=gamvec(i-nnt+1)/4
 647             DU2fric(ind+1)=0.0d0
 648             ind=ind+2
 649           else
 650             DU2fric(ind)=0.0d0
 651             ind=ind+1
 652           endif
 653         enddo
 654       enddo
 655       if (lprn) then
 656       write(iout,*)"The upper part of the five-diagonal friction matrix"
 657       do ichain=1,nchain
 658         write (iout,'(a,i5)') 'Chain',ichain
 659         innt=iposd_chain(ichain)
 660         inct=iposd_chain(ichain)+dimen_chain(ichain)-1
 661         do i=innt,inct
 662           if (i.lt.inct-1) then
 663             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
 664      &       DU2fric(i)
 665           else if (i.eq.inct-1) then
 666             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
 667           else
 668             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
 669           endif
 670         enddo
 671       enddo
 672       endif
 673    10 continue
 674 #else
 675 c The friction matrix
 676       do k=1,dimen
 677        do i=1,dimen
 678          dtdi=0.0d0
 679          do j=1,dimen1
 680            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
 681          enddo
 682          fricmat(k,i)=dtdi
 683        enddo
 684       enddo
 685
 686       if (lprn) then
 687         write (iout,'(//a)') "Matrix fricmat"
 688         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
 689       endif
 690       if (lang.eq.2 .or. lang.eq.3) then
 691 c Mass-scale the friction matrix if non-direct integration will be performed
 692       do i=1,dimen
 693         do j=1,dimen
 694           Ginvfric(i,j)=0.0d0
 695           do k=1,dimen
 696             do l=1,dimen
 697               Ginvfric(i,j)=Ginvfric(i,j)+
 698      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
 699             enddo
 700           enddo
 701         enddo
 702       enddo
 703 c Diagonalize the friction matrix
 704       ind=0
 705       do i=1,dimen
 706         do j=1,i
 707           ind=ind+1
 708           Ghalf(ind)=Ginvfric(i,j)
 709         enddo
 710       enddo
 711       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 712      &  ierr,iwork)
 713       if (lprn) then
 714         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 715      &    " mass-scaled friction matrix"
 716         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 717       endif
 718 c Precompute matrices for tinker stochastic integrator
 719 #ifndef LANG0
 720       do i=1,dimen
 721         do j=1,dimen
 722           mt1(i,j)=0.0d0
 723           mt2(i,j)=0.0d0
 724           do k=1,dimen
 725             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
 726             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
 727           enddo
 728           mt3(j,i)=mt1(i,j)
 729         enddo
 730       enddo
 731 #endif
 732       else if (lang.eq.4) then
 733 c Diagonalize the friction matrix
 734       ind=0
 735       do i=1,dimen
 736         do j=1,i
 737           ind=ind+1
 738           Ghalf(ind)=fricmat(i,j)
 739         enddo
 740       enddo
 741       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 742      &  ierr,iwork)
 743       if (lprn) then
 744         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 745      &    " friction matrix"
 746         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 747       endif
 748 c Determine the number of zero eigenvalues of the friction matrix
 749       nzero=max0(dimen-dimen1,0)
 750 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
 751 c        nzero=nzero+1
 752 c      enddo
 753       write (iout,*) "Number of zero eigenvalues:",nzero
 754       do i=1,dimen
 755         do j=1,dimen
 756           fricmat(i,j)=0.0d0
 757           do k=nzero+1,dimen
 758             fricmat(i,j)=fricmat(i,j)
 759      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k)
 760           enddo
 761         enddo
 762       enddo
 763       if (lprn) then
 764         write (iout,'(//a)') "Generalized inverse of fricmat"
 765         call matout(dimen,dimen,maxres6,maxres6,fricmat)
 766       endif
 767       endif
 768 #ifdef MPI
 769   10  continue
 770       if (nfgtasks.gt.1) then
 771         if (fg_rank.eq.0) then
 772 c The matching BROADCAST for fg processors is called in ERGASTULUM
 773 #ifdef MPI
 774           time00=MPI_Wtime()
 775 #else
 776           time00=tcpu()
 777 #endif
 778           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
 779 #ifdef MPI
 780           time_Bcast=time_Bcast+MPI_Wtime()-time00
 781 #else
 782           time_Bcast=time_Bcast+tcpu()-time00
 783 #endif
 784 c          print *,"Processor",myrank,
 785 c     &       " BROADCAST iorder in SETUP_FRICMAT"
 786         endif
 787 c       licznik=licznik+1
 788 c        write (iout,*) "setup_fricmat licznik",licznik
 789 #ifdef MPI
 790         time00=MPI_Wtime()
 791 #else
 792         time00=tcpu()
 793 #endif
 794 c Scatter the friction matrix
 795         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
 796      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
 797      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
 798 #ifdef TIMING
 799 #ifdef MPI
 800         time_scatter=time_scatter+MPI_Wtime()-time00
 801         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
 802 #else
 803         time_scatter=time_scatter+tcpu()-time00
 804         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
 805 #endif
 806 #endif
 807         do i=1,dimen
 808           do j=1,2*my_ng_count
 809             fricmat(j,i)=fcopy(i,j)
 810           enddo
 811         enddo
 812 c        write (iout,*) "My chunk of fricmat"
 813 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
 814       endif
 815 #endif
 816 #endif
 817       return
 818       end
 819 c-------------------------------------------------------------------------------
 820       subroutine sdarea(gamvec)
 821 c
 822 c Scale the friction coefficients according to solvent accessible surface areas
 823 c Code adapted from TINKER
 824 c AL 9/3/04
 825 c
 826       implicit none
 827       include 'DIMENSIONS'
 828       include 'COMMON.CONTROL'
 829       include 'COMMON.VAR'
 830       include 'COMMON.MD'
 831 #ifndef LANG0
 832       include 'COMMON.LANGEVIN'
 833 #else
 834 #ifdef FIVEDIAG
 835       include 'COMMON.LANGEVIN.lang0.5diag'
 836 #else
 837       include 'COMMON.LANGEVIN.lang0'
 838 #endif
 839 #endif
 840       include 'COMMON.CHAIN'
 841       include 'COMMON.DERIV'
 842       include 'COMMON.GEO'
 843       include 'COMMON.LOCAL'
 844       include 'COMMON.INTERACT'
 845       include 'COMMON.IOUNITS'
 846       include 'COMMON.NAMES'
 847       double precision radius(maxres2),gamvec(maxres2)
 848       double precision twosix
 849       parameter (twosix=1.122462048309372981d0)
 850       logical lprn /.false./
 851       integer i,j,iti,ind
 852       double precision probe,area,ratio
 853 c
 854 c     determine new friction coefficients every few SD steps
 855 c
 856 c     set the atomic radii to estimates of sigma values
 857 c
 858 c      print *,"Entered sdarea"
 859       probe = 0.0d0
 860
 861       do i=1,2*nres
 862         radius(i)=0.0d0
 863       enddo
 864 c  Load peptide group radii
 865       do i=nnt,nct-1
 866         radius(i)=pstok
 867       enddo
 868 c  Load side chain radii
 869       do i=nnt,nct
 870         iti=itype(i)
 871         if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
 872       enddo
 873 c      do i=1,2*nres
 874 c        write (iout,*) "i",i," radius",radius(i)
 875 c      enddo
 876       do i = 1, 2*nres
 877          radius(i) = radius(i) / twosix
 878          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
 879       end do
 880 c
 881 c     scale atomic friction coefficients by accessible area
 882 c
 883       if (lprn) write (iout,*)
 884      &  "Original gammas, surface areas, scaling factors, new gammas, ",
 885      &  "std's of stochastic forces"
 886       ind=0
 887       do i=nnt,nct-1
 888         if (radius(i).gt.0.0d0) then
 889           call surfatom (i,area,radius)
 890           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
 891           if (lprn) write (iout,'(i5,3f10.5,$)')
 892      &      i,gamvec(ind+1),area,ratio
 893           do j=1,3
 894             ind=ind+1
 895             gamvec(ind) = ratio * gamvec(ind)
 896           enddo
 897           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
 898           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
 899         endif
 900       enddo
 901       do i=nnt,nct
 902         if (radius(i+nres).gt.0.0d0) then
 903           call surfatom (i+nres,area,radius)
 904           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
 905           if (lprn) write (iout,'(i5,3f10.5,$)')
 906      &      i,gamvec(ind+1),area,ratio
 907           do j=1,3
 908             ind=ind+1
 909             gamvec(ind) = ratio * gamvec(ind)
 910           enddo
 911           if (itype(i).ne.ntyp1)
 912      &      stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
 913           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
 914         endif
 915       enddo
 916
 917       return
 918       end