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       include 'COMMON.VAR'
 500       include 'COMMON.CHAIN'
 501       include 'COMMON.DERIV'
 502       include 'COMMON.GEO'
 503       include 'COMMON.LOCAL'
 504       include 'COMMON.INTERACT'
 505       include 'COMMON.MD'
 506 #ifdef FIVEDIAG
 507        include 'COMMON.LAGRANGE.5diag'
 508 #else
 509        include 'COMMON.LAGRANGE'
 510 #endif
 511       include 'COMMON.SETUP'
 512       include 'COMMON.TIME1'
 513 c      integer licznik /0/
 514 c      save licznik
 515 #ifndef LANG0
 516       include 'COMMON.LANGEVIN'
 517 #else
 518 #ifdef FIVEDIAG
 519       include 'COMMON.LANGEVIN.lang0.5diag'
 520 #else
 521       include 'COMMON.LANGEVIN.lang0'
 522 #endif
 523 #endif
 524       include 'COMMON.IOUNITS'
 525       integer IERROR
 526       integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
 527       integer ichain,nind
 528       logical lprn /.false./
 529       double precision dtdi,gamvec(MAXRES2)
 530       common /syfek/ gamvec
 531 #ifndef FIVEDIAG
 532       double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2),
 533      & fcopy(maxres2,maxres2)
 534       double precision work(8*maxres2)
 535       integer iwork(maxres2)
 536       common /przechowalnia/ ginvfric,Ghalf,fcopy
 537 #endif
 538
 539 #ifdef MPI
 540       if (fg_rank.ne.king) goto 10
 541 #endif
 542       ind1=0
 543       do i=nnt,nct-1
 544         ind1=ind1+1
 545         gamvec(ind1)=gamp
 546       enddo
 547 c  Load the friction coefficients corresponding to side chains
 548       m=nct-nnt
 549       if (lprn) write (iout,*) "m",m
 550       ind=0
 551 C      gamsc(ntyp1)=1.0d0
 552       do i=nnt,nct
 553         ind=ind+1
 554         ii = ind+m
 555         iti=itype(i)
 556         gamvec(ii)=gamsc(iabs(iti))
 557       enddo
 558       if (surfarea) call sdarea(gamvec)
 559       if (lprn) then
 560         write (iout,*) "Vector gamvec ii",ii
 561         do i=1,ii
 562           write (iout,'(i5,f10.5)') i, gamvec(i)
 563         enddo
 564       endif
 565 #ifdef FIVEDIAG
 566       DMfric=0.0d0
 567       DU1fric=0.0d0
 568       DU2fric=0.0d0
 569       ind=1
 570       do ichain=1,nchain
 571         innt=chain_border(1,ichain)
 572         inct=chain_border(2,ichain)
 573 c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
 574 c DMfric part
 575         DMfric(ind)=gamvec(innt-nnt+1)/4
 576         if (iabs(itype(innt)).eq.10) then
 577           DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
 578           ind=ind+1
 579         else
 580           DMfric(ind+1)=gamvec(m+innt-nnt+1)
 581           ind=ind+2
 582         endif
 583 c        write (iout,*) "DMfric init ind",ind
 584 c DMfric
 585         do i=innt+1,inct-1
 586           DMfric(ind)=gamvec(i-nnt+1)/2
 587           if (iabs(itype(i)).eq.10) then
 588             DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
 589             ind=ind+1
 590           else
 591             DMfric(ind+1)=gamvec(m+i-nnt+1)
 592             ind=ind+2
 593           endif
 594         enddo
 595 c        write (iout,*) "DMfric endloop ind",ind
 596         if (inct.gt.innt) then
 597           DMfric(ind)=gamvec(inct-1-nnt+1)/4
 598           if (iabs(itype(inct)).eq.10) then
 599             DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
 600             ind=ind+1
 601           else
 602             DMfric(ind+1)=gamvec(inct+m-nnt+1)
 603             ind=ind+2
 604           endif
 605         endif
 606 c        write (iout,*) "DMfric end ind",ind
 607       enddo
 608 c DU1fric part
 609       do ichain=1,nchain
 610         ind=iposd_chain(ichain)
 611         innt=chain_border(1,ichain)
 612         inct=chain_border(2,ichain)
 613         do i=innt,inct
 614           if (iabs(itype(i)).ne.10) then
 615             ind=ind+2
 616           else
 617             DU1fric(ind)=gamvec(i-nnt+1)/4
 618             ind=ind+1
 619           endif
 620         enddo
 621       enddo
 622 c DU2fric part
 623       do ichain=1,nchain
 624         ind=iposd_chain(ichain)
 625         innt=chain_border(1,ichain)
 626         inct=chain_border(2,ichain)
 627         do i=innt,inct-1
 628           if (iabs(itype(i)).ne.10) then
 629             DU2fric(ind)=gamvec(i-nnt+1)/4
 630             DU2fric(ind+1)=0.0d0
 631             ind=ind+2
 632           else
 633             DU2fric(ind)=0.0d0
 634             ind=ind+1
 635           endif
 636         enddo
 637       enddo
 638       if (lprn) then
 639       write(iout,*)"The upper part of the five-diagonal friction matrix"
 640       do ichain=1,nchain
 641         write (iout,'(a,i5)') 'Chain',ichain
 642         innt=iposd_chain(ichain)
 643         inct=iposd_chain(ichain)+dimen_chain(ichain)-1
 644         do i=innt,inct
 645           if (i.lt.inct-1) then
 646             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),
 647      &       DU2fric(i)
 648           else if (i.eq.inct-1) then
 649             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
 650           else
 651             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
 652           endif
 653         enddo
 654       enddo
 655       endif
 656    10 continue
 657 #else
 658 c The friction matrix
 659       do k=1,dimen
 660        do i=1,dimen
 661          dtdi=0.0d0
 662          do j=1,dimen1
 663            dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
 664          enddo
 665          fricmat(k,i)=dtdi
 666        enddo
 667       enddo
 668
 669       if (lprn) then
 670         write (iout,'(//a)') "Matrix fricmat"
 671         call matout2(dimen,dimen,maxres2,maxres2,fricmat)
 672       endif
 673       if (lang.eq.2 .or. lang.eq.3) then
 674 c Mass-scale the friction matrix if non-direct integration will be performed
 675       do i=1,dimen
 676         do j=1,dimen
 677           Ginvfric(i,j)=0.0d0
 678           do k=1,dimen
 679             do l=1,dimen
 680               Ginvfric(i,j)=Ginvfric(i,j)+
 681      &          Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
 682             enddo
 683           enddo
 684         enddo
 685       enddo
 686 c Diagonalize the friction matrix
 687       ind=0
 688       do i=1,dimen
 689         do j=1,i
 690           ind=ind+1
 691           Ghalf(ind)=Ginvfric(i,j)
 692         enddo
 693       enddo
 694       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 695      &  ierr,iwork)
 696       if (lprn) then
 697         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 698      &    " mass-scaled friction matrix"
 699         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 700       endif
 701 c Precompute matrices for tinker stochastic integrator
 702 #ifndef LANG0
 703       do i=1,dimen
 704         do j=1,dimen
 705           mt1(i,j)=0.0d0
 706           mt2(i,j)=0.0d0
 707           do k=1,dimen
 708             mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
 709             mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
 710           enddo
 711           mt3(j,i)=mt1(i,j)
 712         enddo
 713       enddo
 714 #endif
 715       else if (lang.eq.4) then
 716 c Diagonalize the friction matrix
 717       ind=0
 718       do i=1,dimen
 719         do j=1,i
 720           ind=ind+1
 721           Ghalf(ind)=fricmat(i,j)
 722         enddo
 723       enddo
 724       call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec,
 725      &  ierr,iwork)
 726       if (lprn) then
 727         write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",
 728      &    " friction matrix"
 729         call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam)
 730       endif
 731 c Determine the number of zero eigenvalues of the friction matrix
 732       nzero=max0(dimen-dimen1,0)
 733 c      do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
 734 c        nzero=nzero+1
 735 c      enddo
 736       write (iout,*) "Number of zero eigenvalues:",nzero
 737       do i=1,dimen
 738         do j=1,dimen
 739           fricmat(i,j)=0.0d0
 740           do k=nzero+1,dimen
 741             fricmat(i,j)=fricmat(i,j)
 742      &        +fricvec(i,k)*fricvec(j,k)/fricgam(k)
 743           enddo
 744         enddo
 745       enddo
 746       if (lprn) then
 747         write (iout,'(//a)') "Generalized inverse of fricmat"
 748         call matout(dimen,dimen,maxres6,maxres6,fricmat)
 749       endif
 750       endif
 751 #ifdef MPI
 752   10  continue
 753       if (nfgtasks.gt.1) then
 754         if (fg_rank.eq.0) then
 755 c The matching BROADCAST for fg processors is called in ERGASTULUM
 756 #ifdef MPI
 757           time00=MPI_Wtime()
 758 #else
 759           time00=tcpu()
 760 #endif
 761           call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
 762 #ifdef MPI
 763           time_Bcast=time_Bcast+MPI_Wtime()-time00
 764 #else
 765           time_Bcast=time_Bcast+tcpu()-time00
 766 #endif
 767 c          print *,"Processor",myrank,
 768 c     &       " BROADCAST iorder in SETUP_FRICMAT"
 769         endif
 770 c       licznik=licznik+1
 771 c        write (iout,*) "setup_fricmat licznik",licznik
 772 #ifdef MPI
 773         time00=MPI_Wtime()
 774 #else
 775         time00=tcpu()
 776 #endif
 777 c Scatter the friction matrix
 778         call MPI_Scatterv(fricmat(1,1),nginv_counts(0),
 779      &    nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),
 780      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
 781 #ifdef TIMING
 782 #ifdef MPI
 783         time_scatter=time_scatter+MPI_Wtime()-time00
 784         time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
 785 #else
 786         time_scatter=time_scatter+tcpu()-time00
 787         time_scatter_fmat=time_scatter_fmat+tcpu()-time00
 788 #endif
 789 #endif
 790         do i=1,dimen
 791           do j=1,2*my_ng_count
 792             fricmat(j,i)=fcopy(i,j)
 793           enddo
 794         enddo
 795 c        write (iout,*) "My chunk of fricmat"
 796 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
 797       endif
 798 #endif
 799 #endif
 800       return
 801       end
 802 c-------------------------------------------------------------------------------
 803       subroutine sdarea(gamvec)
 804 c
 805 c Scale the friction coefficients according to solvent accessible surface areas
 806 c Code adapted from TINKER
 807 c AL 9/3/04
 808 c
 809       implicit none
 810       include 'DIMENSIONS'
 811       include 'COMMON.CONTROL'
 812       include 'COMMON.VAR'
 813       include 'COMMON.MD'
 814 #ifndef LANG0
 815       include 'COMMON.LANGEVIN'
 816 #else
 817 #ifdef FIVEDIAG
 818       include 'COMMON.LANGEVIN.lang0.5diag'
 819 #else
 820       include 'COMMON.LANGEVIN.lang0'
 821 #endif
 822 #endif
 823       include 'COMMON.CHAIN'
 824       include 'COMMON.DERIV'
 825       include 'COMMON.GEO'
 826       include 'COMMON.LOCAL'
 827       include 'COMMON.INTERACT'
 828       include 'COMMON.IOUNITS'
 829       include 'COMMON.NAMES'
 830       double precision radius(maxres2),gamvec(maxres2)
 831       double precision twosix
 832       parameter (twosix=1.122462048309372981d0)
 833       logical lprn /.false./
 834       integer i,j,iti,ind
 835       double precision probe,area,ratio
 836 c
 837 c     determine new friction coefficients every few SD steps
 838 c
 839 c     set the atomic radii to estimates of sigma values
 840 c
 841 c      print *,"Entered sdarea"
 842       probe = 0.0d0
 843
 844       do i=1,2*nres
 845         radius(i)=0.0d0
 846       enddo
 847 c  Load peptide group radii
 848       do i=nnt,nct-1
 849         radius(i)=pstok
 850       enddo
 851 c  Load side chain radii
 852       do i=nnt,nct
 853         iti=itype(i)
 854         if (iti.ne.ntyp1) radius(i+nres)=restok(iti)
 855       enddo
 856 c      do i=1,2*nres
 857 c        write (iout,*) "i",i," radius",radius(i)
 858 c      enddo
 859       do i = 1, 2*nres
 860          radius(i) = radius(i) / twosix
 861          if (radius(i) .ne. 0.0d0)  radius(i) = radius(i) + probe
 862       end do
 863 c
 864 c     scale atomic friction coefficients by accessible area
 865 c
 866       if (lprn) write (iout,*)
 867      &  "Original gammas, surface areas, scaling factors, new gammas, ",
 868      &  "std's of stochastic forces"
 869       ind=0
 870       do i=nnt,nct-1
 871         if (radius(i).gt.0.0d0) then
 872           call surfatom (i,area,radius)
 873           ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
 874           if (lprn) write (iout,'(i5,3f10.5,$)')
 875      &      i,gamvec(ind+1),area,ratio
 876           do j=1,3
 877             ind=ind+1
 878             gamvec(ind) = ratio * gamvec(ind)
 879           enddo
 880           stdforcp(i)=stdfp*dsqrt(gamvec(ind))
 881           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
 882         endif
 883       enddo
 884       do i=nnt,nct
 885         if (radius(i+nres).gt.0.0d0) then
 886           call surfatom (i+nres,area,radius)
 887           ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
 888           if (lprn) write (iout,'(i5,3f10.5,$)')
 889      &      i,gamvec(ind+1),area,ratio
 890           do j=1,3
 891             ind=ind+1
 892             gamvec(ind) = ratio * gamvec(ind)
 893           enddo
 894           if (itype(i).ne.ntyp1)
 895      &      stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
 896           if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
 897         endif
 898       enddo
 899
 900       return
 901       end