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

   1        subroutine lagrangian
   2 c-------------------------------------------------------------------------
   3 c  This subroutine contains the total lagrangain from which the accelerations
   4 c  are obtained.  For numerical gradient checking, the derivetive of the
   5 c  lagrangian in the velocities and coordinates are calculated seperately
   6 c-------------------------------------------------------------------------
   7        implicit none
   8        include 'DIMENSIONS'
   9 #ifdef MPI
  10        include 'mpif.h'
  11        integer time00
  12 #endif
  13        include 'COMMON.VAR'
  14        include 'COMMON.CHAIN'
  15        include 'COMMON.DERIV'
  16        include 'COMMON.GEO'
  17        include 'COMMON.LOCAL'
  18        include 'COMMON.INTERACT'
  19        include 'COMMON.MD'
  20 #ifdef FIVEDIAG
  21        include 'COMMON.LAGRANGE.5diag'
  22 #else
  23        include 'COMMON.LAGRANGE'
  24 #endif
  25 #ifdef LANG0
  26 #ifdef FIVEDIAG
  27       include 'COMMON.LANGEVIN.lang0.5diag'
  28 #else
  29       include 'COMMON.LANGEVIN.lang0'
  30 #endif
  31 #else
  32        include 'COMMON.LANGEVIN'
  33 #endif
  34        include 'COMMON.IOUNITS'
  35        include 'COMMON.CONTROL'
  36        include 'COMMON.MUCA'
  37        include 'COMMON.TIME1'
  38
  39        integer i,j,ind
  40        double precision zapas(MAXRES6),muca_factor
  41        logical lprn /.false./
  42        integer itime
  43        common /cipiszcze/ itime
  44 #ifdef FIVEDIAG
  45        double precision rs(maxres2_chain),xsolv(maxres2_chain),ip4
  46        double precision aaux(3)
  47        integer nind,innt,inct,inct_prev,ichain,n,mark
  48 #ifdef CHECK5DSOL
  49        double precision rscheck(maxres2_chain),rsold(maxres2_chain)
  50 #endif
  51 #endif
  52
  53 #ifdef TIMING
  54        time00=MPI_Wtime()
  55 #endif
  56 #ifdef FIVEDIAG
  57       call grad_transform
  58       d_a=0.0d0
  59       if (lprn) then
  60         write (iout,*) "Potential forces backbone"
  61         do i=1,nres
  62           write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
  63         enddo
  64         write (iout,*) "Potential forces sidechain"
  65         do i=nnt,nct
  66 !          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
  67           write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
  68         enddo
  69       endif
  70       do ichain=1,nchain
  71         n=dimen_chain(ichain)
  72         innt=iposd_chain(ichain)
  73         do j=1,3
  74           ind=1
  75           do i=chain_border(1,ichain),chain_border(2,ichain)
  76             if (itype(i).eq.10)then
  77               rs(ind)=-gcart(j,i)-gxcart(j,i)
  78               ind=ind+1
  79             else
  80               rs(ind)=-gcart(j,i)
  81               rs(ind+1)=-gxcart(j,i)
  82               ind=ind+2
  83             end if
  84           enddo
  85 #ifdef CHECK5DSOL
  86           rsold=rs
  87 #endif
  88           if (lprn) then
  89             write(iout,*)
  90      &      "RHS of the 5-diag equations system, chain",ichain," j",j
  91             do i=1,n
  92               write(iout,'(i5,f10.5)') i,rs(i)
  93             enddo
  94           endif
  95           call FDISYS (n,DM(innt),DU1(innt),DU2(innt),rs,xsolv)
  96           if (lprn) then
  97             write (iout,*) "Solution of the 5-diagonal equations system"
  98             do i=1,n
  99               write (iout,'(i5,f10.5)') i,xsolv(i)
 100             enddo
 101           endif
 102 #ifdef CHECK5DSOL
 103 ! Check the solution
 104           call fivediagmult(n,DMorig(innt),DU1orig(innt),DU2orig(innt),
 105      &      xsolv,rscheck)
 106           do i=1,n
 107             write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),
 108      &       "ratio",rscheck(i)/rsold(i)
 109           enddo
 110 ! end check
 111 #endif
 112 #undef CHECK5DSOL
 113           ind=1
 114           do i=chain_border(1,ichain),chain_border(2,ichain)
 115             if (itype(i).eq.10) then
 116               d_a(j,i)=xsolv(ind)
 117               ind=ind+1
 118             else
 119               d_a(j,i)=xsolv(ind)
 120               d_a(j,i+nres)=xsolv(ind+1)
 121               ind=ind+2
 122             end if
 123           enddo
 124         enddo ! j
 125       enddo ! ichain
 126       if (lprn) then
 127         write (iout,*) "Acceleration in CA and SC oordinates"
 128         do i=1,nres
 129           write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
 130         enddo
 131         do i=nnt,nct
 132           write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
 133         enddo
 134       endif
 135 C Conevert d_a to virtual-bon-vector basis
 136 #define WLOS
 137 #ifdef WLOS
 138 c      write (iout,*) "WLOS"
 139       if (nnt.eq.1) then
 140         d_a(:,0)=d_a(:,1)
 141       endif
 142       do i=1,nres
 143         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 144           do j=1,3
 145             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 146           enddo
 147         else
 148           do j=1,3
 149             d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 150             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 151           enddo
 152         end if
 153       enddo
 154       d_a(:,nres)=0.0d0
 155       d_a(:,nct)=0.0d0
 156       d_a(:,2*nres)=0.0d0
 157 c      d_a(:,0)=d_a(:,1)
 158 c      d_a(:,1)=0.0d0
 159 c      write (iout,*) "Shifting accelerations"
 160       if (nnt.gt.1) then
 161         d_a(:,0)=d_a(:,1)
 162         d_a(:,1)=0.0d0
 163       endif
 164 #define CHUJ
 165 #ifdef CHUJ
 166       do ichain=2,nchain
 167 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
 168 c     &     chain_border1(1,ichain)
 169         d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain))
 170         d_a(:,chain_border1(1,ichain))=0.0d0
 171       enddo
 172 c      write (iout,*) "Adding accelerations"
 173       do ichain=2,nchain
 174 c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
 175 c     &   chain_border(2,ichain-1)
 176         d_a(:,chain_border1(1,ichain)-1)=
 177      &  d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
 178         d_a(:,chain_border(2,ichain-1))=0.0d0
 179       enddo
 180 #endif
 181 #else
 182       inct_prev=0
 183       do j=1,3
 184         aaux(j)=0.0d0
 185       enddo
 186       do ichain=1,nchain
 187         innt=chain_border(1,ichain)
 188         inct=chain_border(2,ichain)
 189         do j=1,3
 190           d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
 191         enddo
 192         inct_prev=inct+1
 193         do i=innt,inct
 194           if (itype(i).ne.10) then
 195             do j=1,3
 196               d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 197             enddo
 198           endif
 199         enddo
 200         do j=1,3
 201           aaux(j)=d_a(j,inct)
 202         enddo
 203         do i=innt,inct
 204           do j=1,3
 205             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 206           enddo
 207         enddo
 208       enddo
 209 #endif
 210       if (lprn) then
 211         write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
 212         do i=0,nres
 213           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 214         enddo
 215         do i=nnt,nct
 216           write (iout,'(i3,3f10.5,3x,3f10.5)')
 217      &     i,(d_a(j,i+nres),j=1,3)
 218         enddo
 219       endif
 220 #else
 221        do j=1,3
 222          zapas(j)=-gcart(j,0)
 223        enddo
 224       ind=3
 225       if (lprn) then
 226         write (iout,*) "Potential forces backbone"
 227       endif
 228       do i=nnt,nct-1
 229         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
 230      &    i,(-gcart(j,i),j=1,3)
 231         do j=1,3
 232           ind=ind+1
 233           zapas(ind)=-gcart(j,i)
 234         enddo
 235       enddo
 236       if (lprn) write (iout,*) "Potential forces sidechain"
 237       do i=nnt,nct
 238         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 239           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
 240      &       i,(-gcart(j,i),j=1,3)
 241           do j=1,3
 242             ind=ind+1
 243             zapas(ind)=-gxcart(j,i)
 244           enddo
 245         endif
 246       enddo
 247
 248       call ginv_mult(zapas,d_a_work)
 249
 250       do j=1,3
 251         d_a(j,0)=d_a_work(j)
 252       enddo
 253       ind=3
 254       do i=nnt,nct-1
 255         do j=1,3
 256           ind=ind+1
 257           d_a(j,i)=d_a_work(ind)
 258         enddo
 259       enddo
 260       do i=nnt,nct
 261         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 262           do j=1,3
 263             ind=ind+1
 264             d_a(j,i+nres)=d_a_work(ind)
 265           enddo
 266         endif
 267       enddo
 268
 269       if(lmuca) then
 270        imtime=imtime+1
 271        if(mucadyn.gt.0) call muca_update(potE)
 272        factor=muca_factor(potE)*t_bath*Rb
 273
 274 cd       print *,'lmuca ',factor,potE
 275        do j=1,3
 276           d_a(j,0)=d_a(j,0)*factor
 277        enddo
 278        do i=nnt,nct-1
 279          do j=1,3
 280           d_a(j,i)=d_a(j,i)*factor
 281          enddo
 282        enddo
 283        do i=nnt,nct
 284          do j=1,3
 285           d_a(j,i+nres)=d_a(j,i+nres)*factor
 286          enddo
 287        enddo
 288
 289       endif
 290
 291       if (lprn) then
 292         write(iout,*) 'acceleration 3D'
 293         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 294         do i=nnt,nct-1
 295          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 296         enddo
 297         do i=nnt,nct
 298          write (iout,'(i3,3f10.5,3x,3f10.5)')
 299      &     i+nres,(d_a(j,i+nres),j=1,3)
 300         enddo
 301       endif
 302 #endif
 303 #ifdef TIMING
 304       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 305 #endif
 306       return
 307       end
 308 c------------------------------------------------------------------
 309       subroutine setup_MD_matrices
 310       implicit none
 311       include 'DIMENSIONS'
 312 #ifdef MPI
 313       include 'mpif.h'
 314       integer ierror,ierr
 315       double precision time00
 316 #endif
 317       include 'COMMON.CONTROL'
 318       include 'COMMON.SETUP'
 319       include 'COMMON.VAR'
 320       include 'COMMON.CHAIN'
 321       include 'COMMON.DERIV'
 322       include 'COMMON.GEO'
 323       include 'COMMON.LOCAL'
 324       include 'COMMON.INTERACT'
 325       include 'COMMON.MD'
 326 #ifdef FIVEDIAG
 327        include 'COMMON.LAGRANGE.5diag'
 328 #else
 329        include 'COMMON.LAGRANGE'
 330 #endif
 331 #ifndef LANG0
 332       include 'COMMON.LANGEVIN'
 333 #else
 334 #ifdef FIVEDIAG
 335       include 'COMMON.LANGEVIN.lang0.5diag'
 336 #else
 337       include 'COMMON.LANGEVIN.lang0'
 338 #endif
 339 #endif
 340       include 'COMMON.IOUNITS'
 341       include 'COMMON.TIME1'
 342       integer i,j,k,m,m1,ind,ind1,ii,iti,ii1,jj
 343       double precision coeff
 344       logical lprn /.false./
 345       logical osob
 346       double precision dtdi,massvec(maxres2)
 347 #ifdef FIVEDIAG
 348       integer ichain,innt,inct,nind,mark,n
 349       double precision ip4
 350 #else
 351       double precision Gcopy(maxres2,maxres2),Ghalf(mmaxres2),
 352      & sqreig(maxres2)
 353       double precision work(8*maxres6)
 354       integer iwork(maxres6)
 355       common /przechowalnia/ Gcopy,Ghalf
 356 #endif
 357 c
 358 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 359 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 360 c
 361 c Determine the number of degrees of freedom (dimen) and the number of
 362 c sites (dimen1)
 363 #ifdef FIVEDIAG
 364       dimen=0
 365       dimen1=0
 366       do ichain=1,nchain
 367         dimen=dimen+chain_length(ichain)
 368         dimen1=dimen1+2*chain_length(ichain)-1
 369         dimenp=dimenp+chain_length(ichain)-1
 370       enddo
 371       write (iout,*) "Number of Calphas",dimen
 372       write (iout,*) "Number of sidechains",nside
 373       write (iout,*) "Number of peptide groups",dimenp
 374       dimen=dimen+nside ! number of centers
 375       dimen3=3*dimen  ! degrees of freedom
 376       write (iout,*) "Number of centers",dimen
 377       write (iout,*) "Degrees of freedom:",dimen3
 378       ip4=ip/4
 379       ind=1
 380       do ichain=1,nchain
 381         iposd_chain(ichain)=ind
 382         innt=chain_border(1,ichain)
 383         inct=chain_border(2,ichain)
 384         DM(ind)=mp/4+ip4
 385         if (iabs(itype(innt)).eq.10) then
 386           DM(ind)=DM(ind)+msc(10)
 387           ind=ind+1
 388           nind=1
 389         else
 390           DM(ind)=DM(ind)+isc(iabs(itype(innt)))
 391           DM(ind+1)=msc(iabs(itype(innt)))+isc(iabs(itype(innt)))
 392           ind=ind+2
 393           nind=2
 394         endif
 395         write (iout,*) "ind",ind," nind",nind
 396         do i=innt+1,inct-1
 397 !        if (iabs(itype(i)).eq.ntyp1) cycle
 398           DM(ind)=2*ip4+mp/2
 399           if (iabs(itype(i)).eq.10) then
 400             if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10)
 401             ind=ind+1
 402             nind=nind+1
 403           else
 404             DM(ind)=DM(ind)+isc(iabs(itype(i)))
 405             DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
 406             ind=ind+2
 407             nind=nind+2
 408           endif
 409           write (iout,*) "i",i," ind",ind," nind",nind
 410         enddo
 411         if (inct.gt.innt) then
 412           DM(ind)=ip4+mp/4
 413           if (iabs(itype(inct)).eq.10) then
 414             DM(ind)=DM(ind)+msc(10)
 415             ind=ind+1
 416             nind=nind+1
 417           else
 418             DM(ind)=DM(ind)+isc(iabs(itype(inct)))
 419             DM(ind+1)=msc(iabs(itype(inct)))+isc(iabs(itype(inct)))
 420             ind=ind+2
 421             nind=nind+2
 422           endif
 423         endif
 424         write (iout,*) "ind",ind," nind",nind
 425         dimen_chain(ichain)=nind
 426       enddo
 427
 428       do ichain=1,nchain
 429         ind=iposd_chain(ichain)
 430         innt=chain_border(1,ichain)
 431         inct=chain_border(2,ichain)
 432         do i=innt,inct
 433           if (iabs(itype(i)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
 434             DU1(ind)=-isc(iabs(itype(i)))
 435             DU1(ind+1)=0.0d0
 436             ind=ind+2
 437           else
 438             DU1(ind)=mp/4-ip4
 439             ind=ind+1
 440           endif
 441         enddo
 442       enddo
 443
 444       do ichain=1,nchain
 445         ind=iposd_chain(ichain)
 446         innt=chain_border(1,ichain)
 447         inct=chain_border(2,ichain)
 448         do i=innt,inct-1
 449 !       if (iabs(itype(i)).eq.ntyp1) cycle
 450 c          write (iout,*) "i",i," itype",itype(i),ntyp1
 451           if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
 452             DU2(ind)=mp/4-ip4
 453             DU2(ind+1)=0.0d0
 454             ind=ind+2
 455           else
 456             DU2(ind)=0.0d0
 457             DU2(ind+1)=0.0d0
 458             ind=ind+1
 459           endif
 460         enddo
 461       enddo
 462       DMorig=DM
 463       DU1orig=DU1
 464       DU2orig=DU2
 465       if (gmatout) then
 466       write (iout,*)"The upper part of the five-diagonal inertia matrix"
 467       endif
 468       do ichain=1,nchain
 469         if (gmatout) write (iout,'(a,i5)') 'Chain',ichain
 470         n=dimen_chain(ichain)
 471         innt=iposd_chain(ichain)
 472         inct=iposd_chain(ichain)+dimen_chain(ichain)-1
 473         if (gmatout) then
 474         do i=innt,inct
 475           if (i.lt.inct-1) then
 476             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 477           else if (i.eq.inct-1) then
 478             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 479           else
 480             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 481           endif
 482         enddo
 483         endif
 484         call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),
 485      &   DU2(innt:inct-1), MARK)
 486
 487         if (mark.eq.-1) then
 488           write(iout,*)
 489      &   "ERROR: the inertia matrix is not positive definite for chain",
 490      &    ichain
 491 #ifdef MPI
 492          call MPI_Finalize(ierr)
 493 #endif
 494          stop
 495         else if (mark.eq.0) then
 496           write (iout,*)
 497      &     "ERROR: the inertia matrix is singular for chain",ichain
 498 #ifdef MPI
 499           call MPI_Finalize(ierr)
 500 #endif
 501         else if (mark.eq.1) then
 502           if (gmatout) then
 503           write (iout,*) "The transformed five-diagonal inertia matrix"
 504           write (iout,'(a,i5)') 'Chain',ichain
 505           do i=innt,inct
 506             if (i.lt.inct-1) then
 507               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 508             else if (i.eq.inct-1) then
 509               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 510             else
 511               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 512             endif
 513           enddo
 514           endif
 515         endif
 516       enddo
 517 ! Diagonalization of the pentadiagonal matrix
 518 #ifdef TIMING
 519       time00=MPI_Wtime()
 520 #endif
 521 #else
 522       dimen=(nct-nnt+1)+nside
 523       dimen1=(nct-nnt)+(nct-nnt+1)
 524       dimen3=dimen*3
 525       write (iout,*) "Degrees_of_freedom",dimen3
 526 #ifdef MPI
 527       if (nfgtasks.gt.1) then
 528       time00=MPI_Wtime()
 529       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 530       time_Bcast=time_Bcast+MPI_Wtime()-time00
 531       call int_bounds(dimen,igmult_start,igmult_end)
 532       igmult_start=igmult_start-1
 533       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 534      &    ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 535       my_ng_count=igmult_end-igmult_start
 536       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 537      &    MPI_INTEGER,FG_COMM,IERROR)
 538       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
 539      & ' absolute rank',myrank,' igmult_start',igmult_start,
 540      & ' igmult_end',igmult_end,' count',my_ng_count
 541       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 542       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 543       call flush(iout)
 544       else
 545 #endif
 546       igmult_start=1
 547       igmult_end=dimen
 548       my_ng_count=dimen
 549 #ifdef MPI
 550       endif
 551 #endif
 552 c      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 553 c  Zeroing out A and fricmat
 554       do i=1,dimen
 555         do j=1,dimen
 556           A(i,j)=0.0D0
 557         enddo
 558       enddo
 559 c  Diagonal elements of the dC part of A and the respective friction coefficients
 560       ind=1
 561       ind1=0
 562       do i=nnt,nct-1
 563         ind=ind+1
 564         ind1=ind1+1
 565         coeff=0.25d0*IP
 566         massvec(ind1)=mp
 567         Gmat(ind,ind)=coeff
 568         A(ind1,ind)=0.5d0
 569       enddo
 570
 571 c  Off-diagonal elements of the dC part of A
 572       k=3
 573       do i=1,nct-nnt
 574         do j=1,i
 575           A(i,j)=1.0d0
 576         enddo
 577       enddo
 578 c  Diagonal elements of the dX part of A and the respective friction coefficients
 579       m=nct-nnt
 580       m1=nct-nnt+1
 581       ind=0
 582       ind1=0
 583       msc(ntyp1)=1.0d0
 584       do i=nnt,nct
 585         ind=ind+1
 586         ii = ind+m
 587         iti=itype(i)
 588         massvec(ii)=msc(iabs(iti))
 589         if (iti.ne.10 .and. iti.ne.ntyp1) then
 590           ind1=ind1+1
 591           ii1= ind1+m1
 592           A(ii,ii1)=1.0d0
 593           Gmat(ii1,ii1)=ISC(iabs(iti))
 594         endif
 595       enddo
 596 c  Off-diagonal elements of the dX part of A
 597       ind=0
 598       k=nct-nnt
 599       do i=nnt,nct
 600         iti=itype(i)
 601         ind=ind+1
 602         do j=nnt,i
 603           ii = ind
 604           jj = j-nnt+1
 605           A(k+ii,jj)=1.0d0
 606         enddo
 607       enddo
 608       if (gmatout) then
 609         write (iout,*)
 610         write (iout,*) "Vector massvec"
 611         do i=1,dimen1
 612           write (iout,*) i,massvec(i)
 613         enddo
 614         write (iout,'(//a)') "A"
 615         call matout(dimen,dimen1,maxres2,maxres2,A)
 616       endif
 617
 618 c Calculate the G matrix (store in Gmat)
 619       do k=1,dimen
 620        do i=1,dimen
 621          dtdi=0.0d0
 622          do j=1,dimen1
 623            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 624          enddo
 625          Gmat(k,i)=Gmat(k,i)+dtdi
 626        enddo
 627       enddo
 628
 629       if (gmatout) then
 630         write (iout,'(//a)') "Gmat"
 631         call matout(dimen,dimen,maxres2,maxres2,Gmat)
 632       endif
 633       do i=1,dimen
 634         do j=1,dimen
 635           Ginv(i,j)=0.0d0
 636           Gcopy(i,j)=Gmat(i,j)
 637         enddo
 638         Ginv(i,i)=1.0d0
 639       enddo
 640 c Invert the G matrix
 641       call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
 642       if (gmatout) then
 643         write (iout,'(//a)') "Ginv"
 644         call matout(dimen,dimen,maxres2,maxres2,Ginv)
 645       endif
 646 #ifdef MPI
 647       if (nfgtasks.gt.1) then
 648         myginv_ng_count=maxres2*my_ng_count
 649         call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 650      &    nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 651         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 652      &    nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 653         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 654         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 655         call flush(iout)
 656 c        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 657 c     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 658 c     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 659 c        call MPI_Barrier(FG_COMM,IERR)
 660         time00=MPI_Wtime()
 661         call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 662      &    nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 663      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 664 #ifdef TIMING
 665         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 666 #endif
 667         do i=1,dimen
 668           do j=1,2*my_ng_count
 669             ginv(j,i)=gcopy(i,j)
 670           enddo
 671         enddo
 672 c        write (iout,*) "Master's chunk of ginv"
 673 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 674       endif
 675 #endif
 676       if (osob) then
 677         write (iout,*) "The G matrix is singular."
 678         stop
 679       endif
 680 c Compute G**(-1/2) and G**(1/2)
 681       ind=0
 682       do i=1,dimen
 683         do j=1,i
 684           ind=ind+1
 685           Ghalf(ind)=Gmat(i,j)
 686         enddo
 687       enddo
 688       call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
 689      &  ierr,iwork)
 690       if (gmatout) then
 691         write (iout,'(//a)')
 692      &   "Eigenvectors and eigenvalues of the G matrix"
 693         call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
 694       endif
 695       do i=1,dimen
 696         sqreig(i)=dsqrt(Geigen(i))
 697       enddo
 698       do i=1,dimen
 699         do j=1,dimen
 700           Gsqrp(i,j)=0.0d0
 701           Gsqrm(i,j)=0.0d0
 702           Gcopy(i,j)=0.0d0
 703           do k=1,dimen
 704             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 705             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 706             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 707           enddo
 708         enddo
 709       enddo
 710       if (lprn) then
 711         write (iout,*) "Comparison of original and restored G"
 712         do i=1,dimen
 713           do j=1,dimen
 714             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
 715      &        Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 716           enddo
 717         enddo
 718       endif
 719 #endif
 720       return
 721       end
 722 c-------------------------------------------------------------------------------
 723       SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
 724       implicit none
 725       include 'DIMENSIONS'
 726       include 'COMMON.IOUNITS'
 727       double precision A(LM2,LM3),B(LM2)
 728       integer nc,nr,lm2,lm3,ka,kb,kc,n,i,j
 729       KA=1
 730       KC=6
 731     1 KB=MIN0(KC,NC)
 732       WRITE(IOUT,600) (I,I=KA,KB)
 733       WRITE(IOUT,601) (B(I),I=KA,KB)
 734       WRITE(IOUT,602)
 735     2 N=0
 736       DO 3  I=1,NR
 737       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 738       N=N+1
 739       IF(N.LT.10) GO TO 3
 740       WRITE(IOUT,602)
 741       N=0
 742     3 CONTINUE
 743     4 IF (KB.EQ.NC) RETURN
 744       KA=KC+1
 745       KC=KC+6
 746       GO TO 1
 747   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 748   601 FORMAT (/5X,10(1PE11.4))
 749   602 FORMAT (2H  )
 750   603 FORMAT (I5,10F11.5)
 751   604 FORMAT (1H1)
 752       END
 753 c-------------------------------------------------------------------------------
 754       SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
 755       implicit none
 756       include 'DIMENSIONS'
 757       include 'COMMON.IOUNITS'
 758       double precision A(LM2,LM3)
 759       integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j
 760       KA=1
 761       KC=6
 762     1 KB=MIN0(KC,NC)
 763       WRITE(IOUT,600) (I,I=KA,KB)
 764       WRITE(IOUT,602)
 765     2 N=0
 766       DO 3  I=1,NR
 767       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 768       N=N+1
 769       IF(N.LT.10) GO TO 3
 770       WRITE(IOUT,602)
 771       N=0
 772     3 CONTINUE
 773     4 IF (KB.EQ.NC) RETURN
 774       KA=KC+1
 775       KC=KC+6
 776       GO TO 1
 777   600 FORMAT (//5x,9I11)
 778   602 FORMAT (2H  )
 779   603 FORMAT (I5,10F11.3)
 780   604 FORMAT (1H1)
 781       END
 782 c-------------------------------------------------------------------------------
 783       SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
 784       implicit none
 785       include 'DIMENSIONS'
 786       include 'COMMON.IOUNITS'
 787       double precision A(LM2,LM3)
 788       integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j
 789       KA=1
 790       KC=21
 791     1 KB=MIN0(KC,NC)
 792       WRITE(IOUT,600) (I,I=KA,KB)
 793       WRITE(IOUT,602)
 794     2 N=0
 795       DO 3  I=1,NR
 796       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 797       N=N+1
 798       IF(N.LT.3) GO TO 3
 799       WRITE(IOUT,602)
 800       N=0
 801     3 CONTINUE
 802     4 IF (KB.EQ.NC) RETURN
 803       KA=KC+1
 804       KC=KC+21
 805       GO TO 1
 806   600 FORMAT (//5x,7(3I5,2x))
 807   602 FORMAT (2H  )
 808   603 FORMAT (I5,7(3F5.1,2x))
 809   604 FORMAT (1H1)
 810       END
 811 c-------------------------------------------------------------------------------
 812       SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
 813       implicit none
 814       include 'DIMENSIONS'
 815       include 'COMMON.IOUNITS'
 816       double precision A(LM2,LM3)
 817       integer nc,nr,lm2,lm3,ka,kb,kc,i,j,n
 818       KA=1
 819       KC=12
 820     1 KB=MIN0(KC,NC)
 821       WRITE(IOUT,600) (I,I=KA,KB)
 822       WRITE(IOUT,602)
 823     2 N=0
 824       DO 3  I=1,NR
 825       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 826       N=N+1
 827       IF(N.LT.3) GO TO 3
 828       WRITE(IOUT,602)
 829       N=0
 830     3 CONTINUE
 831     4 IF (KB.EQ.NC) RETURN
 832       KA=KC+1
 833       KC=KC+12
 834       GO TO 1
 835   600 FORMAT (//5x,4(3I9,2x))
 836   602 FORMAT (2H  )
 837   603 FORMAT (I5,4(3F9.3,2x))
 838   604 FORMAT (1H1)
 839       END
 840 c-----------------------------------------------------------------------------
 841       SUBROUTINE MATOUTR(N,A)
 842 c Prints the lower fragment of a symmetric matix
 843       implicit none
 844       integer n
 845       double precision a(n*(n+1)/2)
 846       integer i,j,k,nlim,jlim,jlim1
 847       CHARACTER*6 LINE6 / '------' /
 848       CHARACTER*12 LINE12 / '------------' /
 849       double precision B(10)
 850       include 'COMMON.IOUNITS'
 851       DO 1 I=1,N,10
 852       NLIM=MIN0(I+9,N)
 853       WRITE (IOUT,1000) (K,K=I,NLIM)
 854       WRITE (IOUT,1020) LINE6,(LINE12,K=I,NLIM)
 855  1000 FORMAT (/7X,10(I5,2X))
 856  1020 FORMAT (A6,10A7)
 857       DO 2 J=I,N
 858       JLIM=MIN0(J,NLIM)
 859       JLIM1=JLIM-I+1
 860       DO 3 K=I,JLIM
 861     3 B(K-I+1)=A(J*(J-1)/2+K)
 862       WRITE (IOUT,1010) J,(B(K),K=1,JLIM1)
 863     2 CONTINUE
 864     1 CONTINUE
 865  1010 FORMAT (I3,3X,10(F7.2))
 866       RETURN
 867       END
 868 #ifdef FIVEDIAG
 869 c---------------------------------------------------------------------------
 870       subroutine fivediagmult(n,DM,DU1,DU2,x,y)
 871       implicit none
 872       integer n
 873       double precision DM(n),DU1(n),DU2(n),x(n),y(n)
 874       integer i
 875       y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
 876       y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
 877       do i=3,n-2
 878         y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i)
 879      &      +DU1(i)*x(i+1)+DU2(i)*x(i+2)
 880       enddo
 881       y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1)
 882      & +DU1(n-1)*x(n)
 883       y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
 884       return
 885       end
 886 c---------------------------------------------------------------------------
 887       subroutine fivediaginv_mult(ndim,forces,d_a_vec)
 888       implicit none
 889       include 'DIMENSIONS'
 890       include 'COMMON.CHAIN'
 891       include 'COMMON.IOUNITS'
 892       include 'COMMON.LAGRANGE.5diag'
 893       include 'COMMON.INTERACT'
 894       include 'COMMON.VAR'
 895       integer ndim
 896       double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
 897      &  xsolv(ndim),d_a_vec(6*nres)
 898       integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
 899       accel=0.0d0
 900       do j=1,3
 901 Compute accelerations in Calpha and SC
 902         do ichain=1,nchain
 903           n=dimen_chain(ichain)
 904           iposc=iposd_chain(ichain)
 905           innt=chain_border(1,ichain)
 906           inct=chain_border(2,ichain)
 907           do i=iposc,iposc+n-1
 908             rs(i-iposc+1)=forces(3*(i-1)+j)
 909           enddo
 910 #ifdef DEBUG
 911           write (iout,*) "j",j," chain",ichain
 912           write (iout,*) "rs"
 913           write (iout,'(f10.5)') (rs(i),i=1,n)
 914 #endif
 915           call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
 916 #ifdef DEBUG
 917           write (iout,*) "xsolv"
 918           write (iout,'(f10.5)') (xsolv(i),i=1,n)
 919 #endif
 920           ind=1
 921           do i=innt,inct
 922             if (itype(i).eq.10)then
 923               accel(j,i)=xsolv(ind)
 924               ind=ind+1
 925             else
 926               accel(j,i)=xsolv(ind)
 927               accel(j,i+nres)=xsolv(ind+1)
 928               ind=ind+2
 929             end if
 930           enddo
 931         enddo
 932       enddo
 933 C Convert d_a to virtual-bon-vector basis
 934 #ifdef DEBUG
 935       write (iout,*) "accel in CA-SC basis"
 936       do i=1,nres
 937         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
 938      &      (accel(j,i+nres),j=1,3)
 939       enddo
 940       write (iout,*) "nnt",nnt
 941 #endif
 942       if (nnt.eq.1) then
 943         accel(:,0)=accel(:,1)
 944       endif
 945       do i=1,nres
 946         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 947           do j=1,3
 948             accel(j,i)=accel(j,i+1)-accel(j,i)
 949           enddo
 950         else
 951           do j=1,3
 952             accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
 953             accel(j,i)=accel(j,i+1)-accel(j,i)
 954           enddo
 955         end if
 956       enddo
 957       accel(:,nres)=0.0d0
 958       accel(:,nct)=0.0d0
 959       accel(:,2*nres)=0.0d0
 960       if (nnt.gt.1) then
 961         accel(:,0)=accel(:,1)
 962         accel(:,1)=0.0d0
 963       endif
 964       do ichain=2,nchain
 965         accel(:,chain_border1(1,ichain)-1)=
 966      &    accel(:,chain_border1(1,ichain))
 967         accel(:,chain_border1(1,ichain))=0.0d0
 968       enddo
 969       do ichain=2,nchain
 970         accel(:,chain_border1(1,ichain)-1)=
 971      &  accel(:,chain_border1(1,ichain)-1)
 972      &   +accel(:,chain_border(2,ichain-1))
 973         accel(:,chain_border(2,ichain-1))=0.0d0
 974       enddo
 975 #ifdef DEBUG
 976       write (iout,*) "accel in fivediaginv_mult: 1"
 977       do i=0,2*nres
 978         write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
 979       enddo
 980 #endif
 981       do j=1,3
 982         d_a_vec(j)=accel(j,0)
 983       enddo
 984       ind=3
 985       do i=nnt,nct-1
 986         do j=1,3
 987           d_a_vec(ind+j)=accel(j,i)
 988         enddo
 989         ind=ind+3
 990       enddo
 991       do i=nnt,nct
 992         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 993           do j=1,3
 994             d_a_vec(ind+j)=accel(j,i+nres)
 995           enddo
 996           ind=ind+3
 997         endif
 998       enddo
 999 #ifdef DEBUG
1000       write (iout,*) "d_a_vec"
1001       write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
1002 #endif
1003       return
1004       end
1005 #else
1006 c---------------------------------------------------------------------------
1007       SUBROUTINE ginv_mult(z,d_a_tmp)
1008       implicit none
1009       include 'DIMENSIONS'
1010 #ifdef MPI
1011       include 'mpif.h'
1012       integer ierr,ierror
1013 #endif
1014       include 'COMMON.SETUP'
1015       include 'COMMON.TIME1'
1016       include 'COMMON.MD'
1017       include 'COMMON.LAGRANGE'
1018       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
1019      &,time01,zcopy(dimen3)
1020       integer i,j,k,ind
1021 #ifdef MPI
1022       if (nfgtasks.gt.1) then
1023         if (fg_rank.eq.0) then
1024 c The matching BROADCAST for fg processors is called in ERGASTULUM
1025           time00=MPI_Wtime()
1026           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
1027           time_Bcast=time_Bcast+MPI_Wtime()-time00
1028 c          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
1029         endif
1030 c        write (2,*) "time00",time00
1031 c        write (2,*) "Before Scatterv"
1032 c        call flush(2)
1033 c        write (2,*) "Whole z (for FG master)"
1034 c        do i=1,dimen
1035 c          write (2,*) i,z(i)
1036 c        enddo
1037 c        call MPI_Barrier(FG_COMM,IERROR)
1038         time00=MPI_Wtime()
1039         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
1040      &    MPI_DOUBLE_PRECISION,
1041      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1042 c        write (2,*) "My chunk of z"
1043         do i=1,3*my_ng_count
1044           z(i)=zcopy(i)
1045 c          write (2,*) i,z(i)
1046         enddo
1047 c        write (2,*) "After SCATTERV"
1048 c        call flush(2)
1049 c        write (2,*) "MPI_Wtime",MPI_Wtime()
1050         time_scatter=time_scatter+MPI_Wtime()-time00
1051 #ifdef TIMING
1052         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
1053 #endif
1054 c        write (2,*) "time_scatter",time_scatter
1055 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
1056 c     &    my_ng_count
1057 c        call flush(2)
1058         time01=MPI_Wtime()
1059         do k=0,2
1060           do i=1,dimen
1061             ind=(i-1)*3+k+1
1062             temp(ind)=0.0d0
1063             do j=1,my_ng_count
1064 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
1065 c     &         Ginv(i,j),z((j-1)*3+k+1),
1066 c     &          Ginv(i,j)*z((j-1)*3+k+1)
1067 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
1068               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
1069             enddo
1070           enddo
1071         enddo
1072         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1073 c        write (2,*) "Before REDUCE"
1074 c        call flush(2)
1075 c        write (2,*) "z before reduce"
1076 c        do i=1,dimen
1077 c          write (2,*) i,temp(i)
1078 c        enddo
1079         time00=MPI_Wtime()
1080         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
1081      &      MPI_SUM,king,FG_COMM,IERR)
1082         time_reduce=time_reduce+MPI_Wtime()-time00
1083 c        write (2,*) "After REDUCE"
1084 c        call flush(2)
1085       else
1086 #endif
1087 #ifdef TIMING
1088         time01=MPI_Wtime()
1089 #endif
1090         do k=0,2
1091           do i=1,dimen
1092             ind=(i-1)*3+k+1
1093             d_a_tmp(ind)=0.0d0
1094             do j=1,dimen
1095 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
1096 c              call flush(2)
1097 c     &         Ginv(i,j),z((j-1)*3+k+1),
1098 c     &          Ginv(i,j)*z((j-1)*3+k+1)
1099               d_a_tmp(ind)=d_a_tmp(ind)
1100      &                         +Ginv(j,i)*z((j-1)*3+k+1)
1101 c              d_a_tmp(ind)=d_a_tmp(ind)
1102 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
1103             enddo
1104           enddo
1105         enddo
1106 #ifdef TIMING
1107         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
1108 #endif
1109 #ifdef MPI
1110       endif
1111 #endif
1112       return
1113       end
1114 c---------------------------------------------------------------------------
1115 #ifdef GINV_MULT
1116       SUBROUTINE ginv_mult_test(z,d_a_tmp)
1117       implicit none
1118       include 'DIMENSIONS'
1119       include 'COMMON.LAGRANGE'
1120       double precision z(dimen),d_a_tmp(dimen)
1121       double precision ztmp(dimen/3),dtmp(dimen/3)
1122
1123 c      do i=1,dimen
1124 c        d_a_tmp(i)=0.0d0
1125 c        do j=1,dimen
1126 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
1127 c        enddo
1128 c      enddo
1129 c
1130 c      return
1131
1132 !ibm* unroll(3)
1133       do k=0,2
1134        do j=1,dimen/3
1135         ztmp(j)=z((j-1)*3+k+1)
1136        enddo
1137
1138        call alignx(16,ztmp(1))
1139        call alignx(16,dtmp(1))
1140        call alignx(16,Ginv(1,1))
1141
1142        do i=1,dimen/3
1143         dtmp(i)=0.0d0
1144         do j=1,dimen/3
1145            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
1146         enddo
1147        enddo
1148        do i=1,dimen/3
1149         ind=(i-1)*3+k+1
1150         d_a_tmp(ind)=dtmp(i)
1151        enddo
1152       enddo
1153       return
1154       end
1155 #endif
1156 c---------------------------------------------------------------------------
1157       SUBROUTINE fricmat_mult(z,d_a_tmp)
1158       include 'DIMENSIONS'
1159 #ifdef MPI
1160       include 'mpif.h'
1161       integer IERROR
1162 #endif
1163       include 'COMMON.MD'
1164       include 'COMMON.LAGRANGE'
1165       include 'COMMON.IOUNITS'
1166       include 'COMMON.SETUP'
1167       include 'COMMON.TIME1'
1168 #ifndef LANG0
1169       include 'COMMON.LANGEVIN'
1170 #else
1171       include 'COMMON.LANGEVIN.lang0'
1172 #endif
1173       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
1174      &,time01,zcopy(dimen3)
1175 #ifdef MPI
1176       if (nfgtasks.gt.1) then
1177         if (fg_rank.eq.0) then
1178 c The matching BROADCAST for fg processors is called in ERGASTULUM
1179           time00=MPI_Wtime()
1180           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
1181           time_Bcast=time_Bcast+MPI_Wtime()-time00
1182 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
1183         endif
1184 c        call MPI_Barrier(FG_COMM,IERROR)
1185         time00=MPI_Wtime()
1186         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
1187      &    MPI_DOUBLE_PRECISION,
1188      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
1189 c        write (2,*) "My chunk of z"
1190         do i=1,3*my_ng_count
1191           z(i)=zcopy(i)
1192 c          write (2,*) i,z(i)
1193         enddo
1194         time_scatter=time_scatter+MPI_Wtime()-time00
1195 #ifdef TIMING
1196         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
1197 #endif
1198         time01=MPI_Wtime()
1199         do k=0,2
1200           do i=1,dimen
1201             ind=(i-1)*3+k+1
1202             temp(ind)=0.0d0
1203             do j=1,my_ng_count
1204               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
1205             enddo
1206           enddo
1207         enddo
1208         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1209 c        write (2,*) "Before REDUCE"
1210 c        write (2,*) "d_a_tmp before reduce"
1211 c        do i=1,dimen3
1212 c          write (2,*) i,temp(i)
1213 c        enddo
1214 c        call flush(2)
1215         time00=MPI_Wtime()
1216         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
1217      &      MPI_SUM,king,FG_COMM,IERR)
1218         time_reduce=time_reduce+MPI_Wtime()-time00
1219 c        write (2,*) "After REDUCE"
1220 c        call flush(2)
1221       else
1222 #endif
1223 #ifdef TIMING
1224         time01=MPI_Wtime()
1225 #endif
1226         do k=0,2
1227          do i=1,dimen
1228           ind=(i-1)*3+k+1
1229           d_a_tmp(ind)=0.0d0
1230           do j=1,dimen
1231              d_a_tmp(ind)=d_a_tmp(ind)
1232      &                           -fricmat(j,i)*z((j-1)*3+k+1)
1233           enddo
1234          enddo
1235         enddo
1236 #ifdef TIMING
1237         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
1238 #endif
1239 #ifdef MPI
1240       endif
1241 #endif
1242 c      write (iout,*) "Vector d_a"
1243 c      do i=1,dimen3
1244 c        write (2,*) i,d_a_tmp(i)
1245 c      enddo
1246       return
1247       end
1248 #endif