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=nnt,nct
  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       do ichain=2,nchain
 165 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
 166 c     &     chain_border1(1,ichain)
 167         d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain))
 168         d_a(:,chain_border1(1,ichain))=0.0d0
 169       enddo
 170 c      write (iout,*) "Adding accelerations"
 171       do ichain=2,nchain
 172 c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
 173 c     &   chain_border(2,ichain-1)
 174         d_a(:,chain_border1(1,ichain)-1)=
 175      &  d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1))
 176         d_a(:,chain_border(2,ichain-1))=0.0d0
 177       enddo
 178 #else
 179       inct_prev=0
 180       do j=1,3
 181         aaux(j)=0.0d0
 182       enddo
 183       do ichain=1,nchain
 184         innt=chain_border(1,ichain)
 185         inct=chain_border(2,ichain)
 186         do j=1,3
 187           d_a(j,inct_prev)=d_a(j,innt)-aaux(j)
 188         enddo
 189         inct_prev=inct+1
 190 #ifdef DEBUG
 191         do i=innt,inct
 192           if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 193             do j=1,3
 194               d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 195             enddo
 196           else
 197             do j=1,3
 198               d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 199               d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 200             enddo
 201           end if
 202         enddo
 203 #else
 204         do i=innt,inct
 205           if (itype(i).ne.10) then
 206             do j=1,3
 207               d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
 208             enddo
 209           endif
 210         enddo
 211         do j=1,3
 212           aaux(j)=d_a(j,inct)
 213         enddo
 214         do i=innt,inct
 215           do j=1,3
 216             d_a(j,i)=d_a(j,i+1)-d_a(j,i)
 217           enddo
 218         enddo
 219 #endif
 220       enddo
 221 #endif
 222       if (lprn) then
 223         write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX'
 224         do i=0,nct-1
 225           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 226         enddo
 227         do i=nnt,nct
 228           write (iout,'(i3,3f10.5,3x,3f10.5)')
 229      &     i+nres,(d_a(j,i+nres),j=1,3)
 230         enddo
 231       endif
 232 #else
 233        do j=1,3
 234          zapas(j)=-gcart(j,0)
 235        enddo
 236       ind=3
 237       if (lprn) then
 238         write (iout,*) "Potential forces backbone"
 239       endif
 240       do i=nnt,nct-1
 241         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
 242      &    i,(-gcart(j,i),j=1,3)
 243         do j=1,3
 244           ind=ind+1
 245           zapas(ind)=-gcart(j,i)
 246         enddo
 247       enddo
 248       if (lprn) write (iout,*) "Potential forces sidechain"
 249       do i=nnt,nct
 250         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 251           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
 252      &       i,(-gcart(j,i),j=1,3)
 253           do j=1,3
 254             ind=ind+1
 255             zapas(ind)=-gxcart(j,i)
 256           enddo
 257         endif
 258       enddo
 259
 260       call ginv_mult(zapas,d_a_work)
 261
 262       do j=1,3
 263         d_a(j,0)=d_a_work(j)
 264       enddo
 265       ind=3
 266       do i=nnt,nct-1
 267         do j=1,3
 268           ind=ind+1
 269           d_a(j,i)=d_a_work(ind)
 270         enddo
 271       enddo
 272       do i=nnt,nct
 273         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 274           do j=1,3
 275             ind=ind+1
 276             d_a(j,i+nres)=d_a_work(ind)
 277           enddo
 278         endif
 279       enddo
 280
 281       if(lmuca) then
 282        imtime=imtime+1
 283        if(mucadyn.gt.0) call muca_update(potE)
 284        factor=muca_factor(potE)*t_bath*Rb
 285
 286 cd       print *,'lmuca ',factor,potE
 287        do j=1,3
 288           d_a(j,0)=d_a(j,0)*factor
 289        enddo
 290        do i=nnt,nct-1
 291          do j=1,3
 292           d_a(j,i)=d_a(j,i)*factor
 293          enddo
 294        enddo
 295        do i=nnt,nct
 296          do j=1,3
 297           d_a(j,i+nres)=d_a(j,i+nres)*factor
 298          enddo
 299        enddo
 300
 301       endif
 302
 303       if (lprn) then
 304         write(iout,*) 'acceleration 3D'
 305         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 306         do i=nnt,nct-1
 307          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 308         enddo
 309         do i=nnt,nct
 310          write (iout,'(i3,3f10.5,3x,3f10.5)')
 311      &     i+nres,(d_a(j,i+nres),j=1,3)
 312         enddo
 313       endif
 314 #endif
 315 #ifdef TIMING
 316       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 317 #endif
 318       return
 319       end
 320 c------------------------------------------------------------------
 321       subroutine setup_MD_matrices
 322       implicit none
 323       include 'DIMENSIONS'
 324 #ifdef MPI
 325       include 'mpif.h'
 326       integer ierror,ierr
 327       double precision time00
 328 #endif
 329       include 'COMMON.CONTROL'
 330       include 'COMMON.SETUP'
 331       include 'COMMON.VAR'
 332       include 'COMMON.CHAIN'
 333       include 'COMMON.DERIV'
 334       include 'COMMON.GEO'
 335       include 'COMMON.LOCAL'
 336       include 'COMMON.INTERACT'
 337       include 'COMMON.MD'
 338 #ifdef FIVEDIAG
 339        include 'COMMON.LAGRANGE.5diag'
 340 #else
 341        include 'COMMON.LAGRANGE'
 342 #endif
 343 #ifndef LANG0
 344       include 'COMMON.LANGEVIN'
 345 #else
 346 #ifdef FIVEDIAG
 347       include 'COMMON.LANGEVIN.lang0.5diag'
 348 #else
 349       include 'COMMON.LANGEVIN.lang0'
 350 #endif
 351 #endif
 352       include 'COMMON.IOUNITS'
 353       include 'COMMON.TIME1'
 354       integer i,j,k,m,m1,ind,ind1,ii,iti,ii1,jj
 355       double precision coeff
 356       logical lprn /.false./
 357       logical osob
 358       double precision dtdi,massvec(maxres2)
 359 #ifdef FIVEDIAG
 360       integer ichain,innt,inct,nind,mark,n
 361       double precision ip4
 362 #else
 363       double precision Gcopy(maxres2,maxres2),Ghalf(mmaxres2),
 364      & sqreig(maxres2)
 365       double precision work(8*maxres6)
 366       integer iwork(maxres6)
 367       common /przechowalnia/ Gcopy,Ghalf
 368 #endif
 369 c
 370 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 371 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 372 c
 373 c Determine the number of degrees of freedom (dimen) and the number of
 374 c sites (dimen1)
 375 #ifdef FIVEDIAG
 376       dimen=0
 377       dimen1=0
 378       do ichain=1,nchain
 379         dimen=dimen+chain_length(ichain)
 380         dimen1=dimen1+2*chain_length(ichain)-1
 381         dimenp=dimenp+chain_length(ichain)-1
 382       enddo
 383       write (iout,*) "Number of Calphas",dimen
 384       write (iout,*) "Number of sidechains",nside
 385       write (iout,*) "Number of peptide groups",dimenp
 386       dimen=dimen+nside ! number of centers
 387       dimen3=3*dimen  ! degrees of freedom
 388       write (iout,*) "Number of centers",dimen
 389       write (iout,*) "Degrees of freedom:",dimen3
 390       ip4=ip/4
 391       ind=1
 392       do ichain=1,nchain
 393         iposd_chain(ichain)=ind
 394         innt=chain_border(1,ichain)
 395         inct=chain_border(2,ichain)
 396         DM(ind)=mp/4+ip4
 397         if (iabs(itype(innt)).eq.10) then
 398           DM(ind)=DM(ind)+msc(10)
 399           ind=ind+1
 400           nind=1
 401         else
 402           DM(ind)=DM(ind)+isc(iabs(itype(innt)))
 403           DM(ind+1)=msc(iabs(itype(innt)))+isc(iabs(itype(innt)))
 404           ind=ind+2
 405           nind=2
 406         endif
 407         write (iout,*) "ind",ind," nind",nind
 408         do i=innt+1,inct-1
 409 !        if (iabs(itype(i)).eq.ntyp1) cycle
 410           DM(ind)=2*ip4+mp/2
 411           if (iabs(itype(i)).eq.10) then
 412             if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10)
 413             ind=ind+1
 414             nind=nind+1
 415           else
 416             DM(ind)=DM(ind)+isc(iabs(itype(i)))
 417             DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
 418             ind=ind+2
 419             nind=nind+2
 420           endif
 421           write (iout,*) "i",i," ind",ind," nind",nind
 422         enddo
 423         if (inct.gt.innt) then
 424           DM(ind)=ip4+mp/4
 425           if (iabs(itype(inct)).eq.10) then
 426             DM(ind)=DM(ind)+msc(10)
 427             ind=ind+1
 428             nind=nind+1
 429           else
 430             DM(ind)=DM(ind)+isc(iabs(itype(inct)))
 431             DM(ind+1)=msc(iabs(itype(inct)))+isc(iabs(itype(inct)))
 432             ind=ind+2
 433             nind=nind+2
 434           endif
 435         endif
 436         write (iout,*) "ind",ind," nind",nind
 437         dimen_chain(ichain)=nind
 438       enddo
 439
 440       do ichain=1,nchain
 441         ind=iposd_chain(ichain)
 442         innt=chain_border(1,ichain)
 443         inct=chain_border(2,ichain)
 444         do i=innt,inct
 445           if (iabs(itype(i)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
 446             DU1(ind)=-isc(iabs(itype(i)))
 447             DU1(ind+1)=0.0d0
 448             ind=ind+2
 449           else
 450             DU1(ind)=mp/4-ip4
 451             ind=ind+1
 452           endif
 453         enddo
 454       enddo
 455
 456       do ichain=1,nchain
 457         ind=iposd_chain(ichain)
 458         innt=chain_border(1,ichain)
 459         inct=chain_border(2,ichain)
 460         do i=innt,inct-1
 461 !       if (iabs(itype(i)).eq.ntyp1) cycle
 462 c          write (iout,*) "i",i," itype",itype(i),ntyp1
 463           if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
 464             DU2(ind)=mp/4-ip4
 465             DU2(ind+1)=0.0d0
 466             ind=ind+2
 467           else
 468             DU2(ind)=0.0d0
 469             DU2(ind+1)=0.0d0
 470             ind=ind+1
 471           endif
 472         enddo
 473       enddo
 474       DMorig=DM
 475       DU1orig=DU1
 476       DU2orig=DU2
 477       if (gmatout) then
 478       write (iout,*)"The upper part of the five-diagonal inertia matrix"
 479       endif
 480       do ichain=1,nchain
 481         if (gmatout) write (iout,'(a,i5)') 'Chain',ichain
 482         n=dimen_chain(ichain)
 483         innt=iposd_chain(ichain)
 484         inct=iposd_chain(ichain)+dimen_chain(ichain)-1
 485         if (gmatout) then
 486         do i=innt,inct
 487           if (i.lt.inct-1) then
 488             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 489           else if (i.eq.inct-1) then
 490             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 491           else
 492             write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 493           endif
 494         enddo
 495         endif
 496         call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),
 497      &   DU2(innt:inct-1), MARK)
 498
 499         if (mark.eq.-1) then
 500           write(iout,*)
 501      &   "ERROR: the inertia matrix is not positive definite for chain",
 502      &    ichain
 503 #ifdef MPI
 504          call MPI_Finalize(ierr)
 505 #endif
 506          stop
 507         else if (mark.eq.0) then
 508           write (iout,*)
 509      &     "ERROR: the inertia matrix is singular for chain",ichain
 510 #ifdef MPI
 511           call MPI_Finalize(ierr)
 512 #endif
 513         else if (mark.eq.1) then
 514           if (gmatout) then
 515           write (iout,*) "The transformed five-diagonal inertia matrix"
 516           write (iout,'(a,i5)') 'Chain',ichain
 517           do i=innt,inct
 518             if (i.lt.inct-1) then
 519               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i)
 520             else if (i.eq.inct-1) then
 521               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i)
 522             else
 523               write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i)
 524             endif
 525           enddo
 526           endif
 527         endif
 528       enddo
 529 ! Diagonalization of the pentadiagonal matrix
 530 #ifdef TIMING
 531       time00=MPI_Wtime()
 532 #endif
 533 #else
 534       dimen=(nct-nnt+1)+nside
 535       dimen1=(nct-nnt)+(nct-nnt+1)
 536       dimen3=dimen*3
 537       write (iout,*) "Degrees_of_freedom",dimen3
 538 #ifdef MPI
 539       if (nfgtasks.gt.1) then
 540       time00=MPI_Wtime()
 541       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 542       time_Bcast=time_Bcast+MPI_Wtime()-time00
 543       call int_bounds(dimen,igmult_start,igmult_end)
 544       igmult_start=igmult_start-1
 545       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 546      &    ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 547       my_ng_count=igmult_end-igmult_start
 548       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 549      &    MPI_INTEGER,FG_COMM,IERROR)
 550       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
 551      & ' absolute rank',myrank,' igmult_start',igmult_start,
 552      & ' igmult_end',igmult_end,' count',my_ng_count
 553       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 554       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 555       call flush(iout)
 556       else
 557 #endif
 558       igmult_start=1
 559       igmult_end=dimen
 560       my_ng_count=dimen
 561 #ifdef MPI
 562       endif
 563 #endif
 564 c      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 565 c  Zeroing out A and fricmat
 566       do i=1,dimen
 567         do j=1,dimen
 568           A(i,j)=0.0D0
 569         enddo
 570       enddo
 571 c  Diagonal elements of the dC part of A and the respective friction coefficients
 572       ind=1
 573       ind1=0
 574       do i=nnt,nct-1
 575         ind=ind+1
 576         ind1=ind1+1
 577         coeff=0.25d0*IP
 578         massvec(ind1)=mp
 579         Gmat(ind,ind)=coeff
 580         A(ind1,ind)=0.5d0
 581       enddo
 582
 583 c  Off-diagonal elements of the dC part of A
 584       k=3
 585       do i=1,nct-nnt
 586         do j=1,i
 587           A(i,j)=1.0d0
 588         enddo
 589       enddo
 590 c  Diagonal elements of the dX part of A and the respective friction coefficients
 591       m=nct-nnt
 592       m1=nct-nnt+1
 593       ind=0
 594       ind1=0
 595       msc(ntyp1)=1.0d0
 596       do i=nnt,nct
 597         ind=ind+1
 598         ii = ind+m
 599         iti=itype(i)
 600         massvec(ii)=msc(iabs(iti))
 601         if (iti.ne.10 .and. iti.ne.ntyp1) then
 602           ind1=ind1+1
 603           ii1= ind1+m1
 604           A(ii,ii1)=1.0d0
 605           Gmat(ii1,ii1)=ISC(iabs(iti))
 606         endif
 607       enddo
 608 c  Off-diagonal elements of the dX part of A
 609       ind=0
 610       k=nct-nnt
 611       do i=nnt,nct
 612         iti=itype(i)
 613         ind=ind+1
 614         do j=nnt,i
 615           ii = ind
 616           jj = j-nnt+1
 617           A(k+ii,jj)=1.0d0
 618         enddo
 619       enddo
 620       if (gmatout) then
 621         write (iout,*)
 622         write (iout,*) "Vector massvec"
 623         do i=1,dimen1
 624           write (iout,*) i,massvec(i)
 625         enddo
 626         write (iout,'(//a)') "A"
 627         call matout(dimen,dimen1,maxres2,maxres2,A)
 628       endif
 629
 630 c Calculate the G matrix (store in Gmat)
 631       do k=1,dimen
 632        do i=1,dimen
 633          dtdi=0.0d0
 634          do j=1,dimen1
 635            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 636          enddo
 637          Gmat(k,i)=Gmat(k,i)+dtdi
 638        enddo
 639       enddo
 640
 641       if (gmatout) then
 642         write (iout,'(//a)') "Gmat"
 643         call matout(dimen,dimen,maxres2,maxres2,Gmat)
 644       endif
 645       do i=1,dimen
 646         do j=1,dimen
 647           Ginv(i,j)=0.0d0
 648           Gcopy(i,j)=Gmat(i,j)
 649         enddo
 650         Ginv(i,i)=1.0d0
 651       enddo
 652 c Invert the G matrix
 653       call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
 654       if (gmatout) then
 655         write (iout,'(//a)') "Ginv"
 656         call matout(dimen,dimen,maxres2,maxres2,Ginv)
 657       endif
 658 #ifdef MPI
 659       if (nfgtasks.gt.1) then
 660         myginv_ng_count=maxres2*my_ng_count
 661         call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 662      &    nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 663         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 664      &    nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 665         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 666         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 667         call flush(iout)
 668 c        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 669 c     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 670 c     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 671 c        call MPI_Barrier(FG_COMM,IERR)
 672         time00=MPI_Wtime()
 673         call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 674      &    nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 675      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 676 #ifdef TIMING
 677         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 678 #endif
 679         do i=1,dimen
 680           do j=1,2*my_ng_count
 681             ginv(j,i)=gcopy(i,j)
 682           enddo
 683         enddo
 684 c        write (iout,*) "Master's chunk of ginv"
 685 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 686       endif
 687 #endif
 688       if (osob) then
 689         write (iout,*) "The G matrix is singular."
 690         stop
 691       endif
 692 c Compute G**(-1/2) and G**(1/2)
 693       ind=0
 694       do i=1,dimen
 695         do j=1,i
 696           ind=ind+1
 697           Ghalf(ind)=Gmat(i,j)
 698         enddo
 699       enddo
 700       call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
 701      &  ierr,iwork)
 702       if (gmatout) then
 703         write (iout,'(//a)')
 704      &   "Eigenvectors and eigenvalues of the G matrix"
 705         call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
 706       endif
 707       do i=1,dimen
 708         sqreig(i)=dsqrt(Geigen(i))
 709       enddo
 710       do i=1,dimen
 711         do j=1,dimen
 712           Gsqrp(i,j)=0.0d0
 713           Gsqrm(i,j)=0.0d0
 714           Gcopy(i,j)=0.0d0
 715           do k=1,dimen
 716             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 717             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 718             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 719           enddo
 720         enddo
 721       enddo
 722       if (lprn) then
 723         write (iout,*) "Comparison of original and restored G"
 724         do i=1,dimen
 725           do j=1,dimen
 726             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
 727      &        Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 728           enddo
 729         enddo
 730       endif
 731 #endif
 732       return
 733       end
 734 c-------------------------------------------------------------------------------
 735       SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
 736       implicit none
 737       include 'DIMENSIONS'
 738       include 'COMMON.IOUNITS'
 739       double precision A(LM2,LM3),B(LM2)
 740       integer nc,nr,lm2,lm3,ka,kb,kc,n,i,j
 741       KA=1
 742       KC=6
 743     1 KB=MIN0(KC,NC)
 744       WRITE(IOUT,600) (I,I=KA,KB)
 745       WRITE(IOUT,601) (B(I),I=KA,KB)
 746       WRITE(IOUT,602)
 747     2 N=0
 748       DO 3  I=1,NR
 749       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 750       N=N+1
 751       IF(N.LT.10) GO TO 3
 752       WRITE(IOUT,602)
 753       N=0
 754     3 CONTINUE
 755     4 IF (KB.EQ.NC) RETURN
 756       KA=KC+1
 757       KC=KC+6
 758       GO TO 1
 759   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 760   601 FORMAT (/5X,10(1PE11.4))
 761   602 FORMAT (2H  )
 762   603 FORMAT (I5,10F11.5)
 763   604 FORMAT (1H1)
 764       END
 765 c-------------------------------------------------------------------------------
 766       SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
 767       implicit none
 768       include 'DIMENSIONS'
 769       include 'COMMON.IOUNITS'
 770       double precision A(LM2,LM3)
 771       integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j
 772       KA=1
 773       KC=6
 774     1 KB=MIN0(KC,NC)
 775       WRITE(IOUT,600) (I,I=KA,KB)
 776       WRITE(IOUT,602)
 777     2 N=0
 778       DO 3  I=1,NR
 779       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 780       N=N+1
 781       IF(N.LT.10) GO TO 3
 782       WRITE(IOUT,602)
 783       N=0
 784     3 CONTINUE
 785     4 IF (KB.EQ.NC) RETURN
 786       KA=KC+1
 787       KC=KC+6
 788       GO TO 1
 789   600 FORMAT (//5x,9I11)
 790   602 FORMAT (2H  )
 791   603 FORMAT (I5,10F11.3)
 792   604 FORMAT (1H1)
 793       END
 794 c-------------------------------------------------------------------------------
 795       SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
 796       implicit none
 797       include 'DIMENSIONS'
 798       include 'COMMON.IOUNITS'
 799       double precision A(LM2,LM3)
 800       integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j
 801       KA=1
 802       KC=21
 803     1 KB=MIN0(KC,NC)
 804       WRITE(IOUT,600) (I,I=KA,KB)
 805       WRITE(IOUT,602)
 806     2 N=0
 807       DO 3  I=1,NR
 808       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 809       N=N+1
 810       IF(N.LT.3) GO TO 3
 811       WRITE(IOUT,602)
 812       N=0
 813     3 CONTINUE
 814     4 IF (KB.EQ.NC) RETURN
 815       KA=KC+1
 816       KC=KC+21
 817       GO TO 1
 818   600 FORMAT (//5x,7(3I5,2x))
 819   602 FORMAT (2H  )
 820   603 FORMAT (I5,7(3F5.1,2x))
 821   604 FORMAT (1H1)
 822       END
 823 c-------------------------------------------------------------------------------
 824       SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
 825       implicit none
 826       include 'DIMENSIONS'
 827       include 'COMMON.IOUNITS'
 828       double precision A(LM2,LM3)
 829       integer nc,nr,lm2,lm3,ka,kb,kc,i,j,n
 830       KA=1
 831       KC=12
 832     1 KB=MIN0(KC,NC)
 833       WRITE(IOUT,600) (I,I=KA,KB)
 834       WRITE(IOUT,602)
 835     2 N=0
 836       DO 3  I=1,NR
 837       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 838       N=N+1
 839       IF(N.LT.3) GO TO 3
 840       WRITE(IOUT,602)
 841       N=0
 842     3 CONTINUE
 843     4 IF (KB.EQ.NC) RETURN
 844       KA=KC+1
 845       KC=KC+12
 846       GO TO 1
 847   600 FORMAT (//5x,4(3I9,2x))
 848   602 FORMAT (2H  )
 849   603 FORMAT (I5,4(3F9.3,2x))
 850   604 FORMAT (1H1)
 851       END
 852 c-----------------------------------------------------------------------------
 853       SUBROUTINE MATOUTR(N,A)
 854 c Prints the lower fragment of a symmetric matix
 855       implicit none
 856       integer n
 857       double precision a(n*(n+1)/2)
 858       integer i,j,k,nlim,jlim,jlim1
 859       CHARACTER*6 LINE6 / '------' /
 860       CHARACTER*12 LINE12 / '------------' /
 861       double precision B(10)
 862       include 'COMMON.IOUNITS'
 863       DO 1 I=1,N,10
 864       NLIM=MIN0(I+9,N)
 865       WRITE (IOUT,1000) (K,K=I,NLIM)
 866       WRITE (IOUT,1020) LINE6,(LINE12,K=I,NLIM)
 867  1000 FORMAT (/7X,10(I5,2X))
 868  1020 FORMAT (A6,10A7)
 869       DO 2 J=I,N
 870       JLIM=MIN0(J,NLIM)
 871       JLIM1=JLIM-I+1
 872       DO 3 K=I,JLIM
 873     3 B(K-I+1)=A(J*(J-1)/2+K)
 874       WRITE (IOUT,1010) J,(B(K),K=1,JLIM1)
 875     2 CONTINUE
 876     1 CONTINUE
 877  1010 FORMAT (I3,3X,10(F7.2))
 878       RETURN
 879       END
 880 #ifdef FIVEDIAG
 881 c---------------------------------------------------------------------------
 882       subroutine fivediagmult(n,DM,DU1,DU2,x,y)
 883       implicit none
 884       integer n
 885       double precision DM(n),DU1(n),DU2(n),x(n),y(n)
 886       integer i
 887       y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
 888       y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
 889       do i=3,n-2
 890         y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i)
 891      &      +DU1(i)*x(i+1)+DU2(i)*x(i+2)
 892       enddo
 893       y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1)
 894      & +DU1(n-1)*x(n)
 895       y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
 896       return
 897       end
 898 c---------------------------------------------------------------------------
 899       subroutine fivediaginv_mult(ndim,forces,d_a_vec)
 900       implicit none
 901       include 'DIMENSIONS'
 902       include 'COMMON.CHAIN'
 903       include 'COMMON.IOUNITS'
 904       include 'COMMON.LAGRANGE.5diag'
 905       include 'COMMON.INTERACT'
 906       integer ndim
 907       double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
 908      &  xsolv(ndim),d_a_vec(6*nres)
 909       integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
 910       do j=1,3
 911 Compute accelerations in Calpha and SC
 912         do ichain=1,nchain
 913           n=dimen_chain(ichain)
 914           iposc=iposd_chain(ichain)
 915           innt=chain_border(1,ichain)
 916           inct=chain_border(2,ichain)
 917           do i=iposc,iposc+n-1
 918             rs(i)=forces(3*(i-1)+j)
 919           enddo
 920           call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
 921           ind=1
 922           do i=innt,inct
 923             if (itype(i).eq.10)then
 924               accel(j,i)=xsolv(ind)
 925               ind=ind+1
 926             else
 927               accel(j,i)=xsolv(ind)
 928               accel(j,i+nres)=xsolv(ind+1)
 929               ind=ind+2
 930             end if
 931           enddo
 932         enddo
 933       enddo
 934 C Conevert d_a to virtual-bon-vector basis
 935 #ifdef DEBUG
 936       write (iout,*) "accel in CA-SC basis"
 937       do i=1,nres
 938         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
 939      &      (accel(j,i+nres),j=1,3)
 940       enddo
 941       write (iout,*) "nnt",nnt
 942 #endif
 943       if (nnt.eq.1) then
 944         accel(:,0)=accel(:,1)
 945       endif
 946       do i=1,nres
 947         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
 948           do j=1,3
 949             accel(j,i)=accel(j,i+1)-accel(j,i)
 950           enddo
 951         else
 952           do j=1,3
 953             accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
 954             accel(j,i)=accel(j,i+1)-accel(j,i)
 955           enddo
 956         end if
 957       enddo
 958       accel(:,nres)=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,dimen3)
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