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