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