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

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