source/unres/src_MD-M/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 real*8 (a-h,o-z)
   8        include 'DIMENSIONS'
   9 #ifdef MPI
  10        include 'mpif.h'
  11 #endif
  12        include 'COMMON.VAR'
  13        include 'COMMON.CHAIN'
  14        include 'COMMON.DERIV'
  15        include 'COMMON.GEO'
  16        include 'COMMON.LOCAL'
  17        include 'COMMON.INTERACT'
  18        include 'COMMON.MD'
  19        include 'COMMON.IOUNITS'
  20        include 'COMMON.CONTROL'
  21        include 'COMMON.MUCA'
  22        include 'COMMON.TIME1'
  23
  24        integer i,j,ind
  25        double precision zapas(MAXRES6),muca_factor
  26        logical lprn /.false./
  27        common /cipiszcze/ itime
  28
  29 #ifdef TIMING
  30        time00=MPI_Wtime()
  31 #endif
  32        do j=1,3
  33          zapas(j)=-gcart(j,0)
  34        enddo
  35       ind=3
  36       if (lprn) then
  37         write (iout,*) "Potential forces backbone"
  38       endif
  39       do i=nnt,nct-1
  40         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
  41      &    i,(-gcart(j,i),j=1,3)
  42         do j=1,3
  43           ind=ind+1
  44           zapas(ind)=-gcart(j,i)
  45         enddo
  46       enddo
  47       if (lprn) write (iout,*) "Potential forces sidechain"
  48       do i=nnt,nct
  49         if (itype(i).ne.10 .and. itype(i).ne.21) then
  50           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)')
  51      &       i,(-gcart(j,i),j=1,3)
  52           do j=1,3
  53             ind=ind+1
  54             zapas(ind)=-gxcart(j,i)
  55           enddo
  56         endif
  57       enddo
  58
  59       call ginv_mult(zapas,d_a_work)
  60
  61       do j=1,3
  62         d_a(j,0)=d_a_work(j)
  63       enddo
  64       ind=3
  65       do i=nnt,nct-1
  66         do j=1,3
  67           ind=ind+1
  68           d_a(j,i)=d_a_work(ind)
  69         enddo
  70       enddo
  71       do i=nnt,nct
  72         if (itype(i).ne.10 .and. itype(i).ne.21) then
  73           do j=1,3
  74             ind=ind+1
  75             d_a(j,i+nres)=d_a_work(ind)
  76           enddo
  77         endif
  78       enddo
  79
  80       if(lmuca) then
  81        imtime=imtime+1
  82        if(mucadyn.gt.0) call muca_update(potE)
  83        factor=muca_factor(potE)*t_bath*Rb
  84
  85 cd       print *,'lmuca ',factor,potE
  86        do j=1,3
  87           d_a(j,0)=d_a(j,0)*factor
  88        enddo
  89        do i=nnt,nct-1
  90          do j=1,3
  91           d_a(j,i)=d_a(j,i)*factor
  92          enddo
  93        enddo
  94        do i=nnt,nct
  95          do j=1,3
  96           d_a(j,i+nres)=d_a(j,i+nres)*factor
  97          enddo
  98        enddo
  99
 100       endif
 101
 102       if (lprn) then
 103         write(iout,*) 'acceleration 3D'
 104         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 105         do i=nnt,nct-1
 106          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 107         enddo
 108         do i=nnt,nct
 109          write (iout,'(i3,3f10.5,3x,3f10.5)')
 110      &     i+nres,(d_a(j,i+nres),j=1,3)
 111         enddo
 112       endif
 113 #ifdef TIMING
 114       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 115 #endif
 116       return
 117       end
 118 c------------------------------------------------------------------
 119       subroutine setup_MD_matrices
 120       implicit real*8 (a-h,o-z)
 121       include 'DIMENSIONS'
 122 #ifdef MPI
 123       include 'mpif.h'
 124       integer ierror
 125 #endif
 126       include 'COMMON.SETUP'
 127       include 'COMMON.VAR'
 128       include 'COMMON.CHAIN'
 129       include 'COMMON.DERIV'
 130       include 'COMMON.GEO'
 131       include 'COMMON.LOCAL'
 132       include 'COMMON.INTERACT'
 133       include 'COMMON.MD'
 134 #ifndef LANG0
 135       include 'COMMON.LANGEVIN'
 136 #else
 137       include 'COMMON.LANGEVIN.lang0'
 138 #endif
 139       include 'COMMON.IOUNITS'
 140       include 'COMMON.TIME1'
 141       integer i,j
 142       logical lprn /.false./
 143       logical osob
 144       double precision dtdi,massvec(maxres2),Gcopy(maxres2,maxres2),
 145      &  Ghalf(mmaxres2),sqreig(maxres2)
 146       double precision work(8*maxres6)
 147       integer iwork(maxres6)
 148       common /przechowalnia/ Gcopy,Ghalf
 149 c
 150 c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 151 c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 152 c
 153 c Determine the number of degrees of freedom (dimen) and the number of
 154 c sites (dimen1)
 155       dimen=(nct-nnt+1)+nside
 156       dimen1=(nct-nnt)+(nct-nnt+1)
 157       dimen3=dimen*3
 158 #ifdef MPI
 159       if (nfgtasks.gt.1) then
 160       time00=MPI_Wtime()
 161       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 162       time_Bcast=time_Bcast+MPI_Wtime()-time00
 163       call int_bounds(dimen,igmult_start,igmult_end)
 164       igmult_start=igmult_start-1
 165       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
 166      &    ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 167       my_ng_count=igmult_end-igmult_start
 168       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
 169      &    MPI_INTEGER,FG_COMM,IERROR)
 170       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,
 171      & ' absolute rank',myrank,' igmult_start',igmult_start,
 172      & ' igmult_end',igmult_end,' count',my_ng_count
 173       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 174       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 175       call flush(iout)
 176       else
 177 #endif
 178       igmult_start=1
 179       igmult_end=dimen
 180       my_ng_count=dimen
 181 #ifdef MPI
 182       endif
 183 #endif
 184 c      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 185 c  Zeroing out A and fricmat
 186       do i=1,dimen
 187         do j=1,dimen
 188           A(i,j)=0.0D0
 189         enddo
 190       enddo
 191 c  Diagonal elements of the dC part of A and the respective friction coefficients
 192       ind=1
 193       ind1=0
 194       do i=nnt,nct-1
 195         ind=ind+1
 196         ind1=ind1+1
 197         coeff=0.25d0*IP
 198         massvec(ind1)=mp
 199         Gmat(ind,ind)=coeff
 200         A(ind1,ind)=0.5d0
 201       enddo
 202
 203 c  Off-diagonal elements of the dC part of A
 204       k=3
 205       do i=1,nct-nnt
 206         do j=1,i
 207           A(i,j)=1.0d0
 208         enddo
 209       enddo
 210 c  Diagonal elements of the dX part of A and the respective friction coefficients
 211       m=nct-nnt
 212       m1=nct-nnt+1
 213       ind=0
 214       ind1=0
 215       msc(21)=1.0d0
 216       do i=nnt,nct
 217         ind=ind+1
 218         ii = ind+m
 219         iti=itype(i)
 220         massvec(ii)=msc(iabs(iti))
 221         if (iti.ne.10 .and. iti.ne.ntyp1) then
 222           ind1=ind1+1
 223           ii1= ind1+m1
 224           A(ii,ii1)=1.0d0
 225           Gmat(ii1,ii1)=ISC(iabs(iti))
 226         endif
 227       enddo
 228 c  Off-diagonal elements of the dX part of A
 229       ind=0
 230       k=nct-nnt
 231       do i=nnt,nct
 232         iti=itype(i)
 233         ind=ind+1
 234         do j=nnt,i
 235           ii = ind
 236           jj = j-nnt+1
 237           A(k+ii,jj)=1.0d0
 238         enddo
 239       enddo
 240       if (lprn) then
 241         write (iout,*)
 242         write (iout,*) "Vector massvec"
 243         do i=1,dimen1
 244           write (iout,*) i,massvec(i)
 245         enddo
 246         write (iout,'(//a)') "A"
 247         call matout(dimen,dimen1,maxres2,maxres2,A)
 248       endif
 249
 250 c Calculate the G matrix (store in Gmat)
 251       do k=1,dimen
 252        do i=1,dimen
 253          dtdi=0.0d0
 254          do j=1,dimen1
 255            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 256          enddo
 257          Gmat(k,i)=Gmat(k,i)+dtdi
 258        enddo
 259       enddo
 260
 261       if (lprn) then
 262         write (iout,'(//a)') "Gmat"
 263         call matout(dimen,dimen,maxres2,maxres2,Gmat)
 264       endif
 265       do i=1,dimen
 266         do j=1,dimen
 267           Ginv(i,j)=0.0d0
 268           Gcopy(i,j)=Gmat(i,j)
 269         enddo
 270         Ginv(i,i)=1.0d0
 271       enddo
 272 c Invert the G matrix
 273       call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
 274       if (lprn) then
 275         write (iout,'(//a)') "Ginv"
 276         call matout(dimen,dimen,maxres2,maxres2,Ginv)
 277       endif
 278 #ifdef MPI
 279       if (nfgtasks.gt.1) then
 280         myginv_ng_count=maxres2*my_ng_count
 281         call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 282      &    nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 283         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 284      &    nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 285         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 286         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 287         call flush(iout)
 288 c        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 289 c     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 290 c     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 291 c        call MPI_Barrier(FG_COMM,IERR)
 292         time00=MPI_Wtime()
 293         call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 294      &    nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 295      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 296 #ifdef TIMING
 297         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 298 #endif
 299         do i=1,dimen
 300           do j=1,2*my_ng_count
 301             ginv(j,i)=gcopy(i,j)
 302           enddo
 303         enddo
 304 c        write (iout,*) "Master's chunk of ginv"
 305 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 306       endif
 307 #endif
 308       if (osob) then
 309         write (iout,*) "The G matrix is singular."
 310         stop
 311       endif
 312 c Compute G**(-1/2) and G**(1/2)
 313       ind=0
 314       do i=1,dimen
 315         do j=1,i
 316           ind=ind+1
 317           Ghalf(ind)=Gmat(i,j)
 318         enddo
 319       enddo
 320       call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
 321      &  ierr,iwork)
 322       if (lprn) then
 323         write (iout,'(//a)')
 324      &   "Eigenvectors and eigenvalues of the G matrix"
 325         call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
 326       endif
 327       do i=1,dimen
 328         sqreig(i)=dsqrt(Geigen(i))
 329       enddo
 330       do i=1,dimen
 331         do j=1,dimen
 332           Gsqrp(i,j)=0.0d0
 333           Gsqrm(i,j)=0.0d0
 334           Gcopy(i,j)=0.0d0
 335           do k=1,dimen
 336             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 337             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 338             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 339           enddo
 340         enddo
 341       enddo
 342       if (lprn) then
 343         write (iout,*) "Comparison of original and restored G"
 344         do i=1,dimen
 345           do j=1,dimen
 346             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
 347      &        Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 348           enddo
 349         enddo
 350       endif
 351       return
 352       end
 353 c-------------------------------------------------------------------------------
 354       SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
 355       implicit real*8 (a-h,o-z)
 356       include 'DIMENSIONS'
 357       include 'COMMON.IOUNITS'
 358       double precision A(LM2,LM3),B(LM2)
 359       KA=1
 360       KC=6
 361     1 KB=MIN0(KC,NC)
 362       WRITE(IOUT,600) (I,I=KA,KB)
 363       WRITE(IOUT,601) (B(I),I=KA,KB)
 364       WRITE(IOUT,602)
 365     2 N=0
 366       DO 3  I=1,NR
 367       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 368       N=N+1
 369       IF(N.LT.10) GO TO 3
 370       WRITE(IOUT,602)
 371       N=0
 372     3 CONTINUE
 373     4 IF (KB.EQ.NC) RETURN
 374       KA=KC+1
 375       KC=KC+6
 376       GO TO 1
 377   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 378   601 FORMAT (/5X,10(1PE11.4))
 379   602 FORMAT (2H  )
 380   603 FORMAT (I5,10F11.5)
 381   604 FORMAT (1H1)
 382       END
 383 c-------------------------------------------------------------------------------
 384       SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
 385       implicit real*8 (a-h,o-z)
 386       include 'DIMENSIONS'
 387       include 'COMMON.IOUNITS'
 388       double precision A(LM2,LM3)
 389       KA=1
 390       KC=6
 391     1 KB=MIN0(KC,NC)
 392       WRITE(IOUT,600) (I,I=KA,KB)
 393       WRITE(IOUT,602)
 394     2 N=0
 395       DO 3  I=1,NR
 396       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 397       N=N+1
 398       IF(N.LT.10) GO TO 3
 399       WRITE(IOUT,602)
 400       N=0
 401     3 CONTINUE
 402     4 IF (KB.EQ.NC) RETURN
 403       KA=KC+1
 404       KC=KC+6
 405       GO TO 1
 406   600 FORMAT (//5x,9I11)
 407   602 FORMAT (2H  )
 408   603 FORMAT (I5,10F11.3)
 409   604 FORMAT (1H1)
 410       END
 411 c-------------------------------------------------------------------------------
 412       SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
 413       implicit real*8 (a-h,o-z)
 414       include 'DIMENSIONS'
 415       include 'COMMON.IOUNITS'
 416       double precision A(LM2,LM3)
 417       KA=1
 418       KC=21
 419     1 KB=MIN0(KC,NC)
 420       WRITE(IOUT,600) (I,I=KA,KB)
 421       WRITE(IOUT,602)
 422     2 N=0
 423       DO 3  I=1,NR
 424       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 425       N=N+1
 426       IF(N.LT.3) GO TO 3
 427       WRITE(IOUT,602)
 428       N=0
 429     3 CONTINUE
 430     4 IF (KB.EQ.NC) RETURN
 431       KA=KC+1
 432       KC=KC+21
 433       GO TO 1
 434   600 FORMAT (//5x,7(3I5,2x))
 435   602 FORMAT (2H  )
 436   603 FORMAT (I5,7(3F5.1,2x))
 437   604 FORMAT (1H1)
 438       END
 439 c-------------------------------------------------------------------------------
 440       SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
 441       implicit real*8 (a-h,o-z)
 442       include 'DIMENSIONS'
 443       include 'COMMON.IOUNITS'
 444       double precision A(LM2,LM3)
 445       KA=1
 446       KC=12
 447     1 KB=MIN0(KC,NC)
 448       WRITE(IOUT,600) (I,I=KA,KB)
 449       WRITE(IOUT,602)
 450     2 N=0
 451       DO 3  I=1,NR
 452       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 453       N=N+1
 454       IF(N.LT.3) GO TO 3
 455       WRITE(IOUT,602)
 456       N=0
 457     3 CONTINUE
 458     4 IF (KB.EQ.NC) RETURN
 459       KA=KC+1
 460       KC=KC+12
 461       GO TO 1
 462   600 FORMAT (//5x,4(3I9,2x))
 463   602 FORMAT (2H  )
 464   603 FORMAT (I5,4(3F9.3,2x))
 465   604 FORMAT (1H1)
 466       END
 467 c---------------------------------------------------------------------------
 468       SUBROUTINE ginv_mult(z,d_a_tmp)
 469       implicit real*8 (a-h,o-z)
 470       include 'DIMENSIONS'
 471 #ifdef MPI
 472       include 'mpif.h'
 473       integer ierr
 474 #endif
 475       include 'COMMON.SETUP'
 476       include 'COMMON.TIME1'
 477       include 'COMMON.MD'
 478       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 479      &time01
 480 #ifdef MPI
 481       if (nfgtasks.gt.1) then
 482         if (fg_rank.eq.0) then
 483 c The matching BROADCAST for fg processors is called in ERGASTULUM
 484           time00=MPI_Wtime()
 485           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
 486           time_Bcast=time_Bcast+MPI_Wtime()-time00
 487 c          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
 488         endif
 489 c        write (2,*) "time00",time00
 490 c        write (2,*) "Before Scatterv"
 491 c        call flush(2)
 492 c        write (2,*) "Whole z (for FG master)"
 493 c        do i=1,dimen
 494 c          write (2,*) i,z(i)
 495 c        enddo
 496 c        call MPI_Barrier(FG_COMM,IERROR)
 497         time00=MPI_Wtime()
 498         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 499      &    MPI_DOUBLE_PRECISION,
 500      &    z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 501 c        write (2,*) "My chunk of z"
 502 c        do i=1,3*my_ng_count
 503 c          write (2,*) i,z(i)
 504 c        enddo
 505 c        write (2,*) "After SCATTERV"
 506 c        call flush(2)
 507 c        write (2,*) "MPI_Wtime",MPI_Wtime()
 508         time_scatter=time_scatter+MPI_Wtime()-time00
 509 #ifdef TIMING
 510         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 511 #endif
 512 c        write (2,*) "time_scatter",time_scatter
 513 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 514 c     &    my_ng_count
 515 c        call flush(2)
 516         time01=MPI_Wtime()
 517         do k=0,2
 518           do i=1,dimen
 519             ind=(i-1)*3+k+1
 520             temp(ind)=0.0d0
 521             do j=1,my_ng_count
 522 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 523 c     &         Ginv(i,j),z((j-1)*3+k+1),
 524 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 525 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 526               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
 527             enddo
 528           enddo
 529         enddo
 530         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 531 c        write (2,*) "Before REDUCE"
 532 c        call flush(2)
 533 c        write (2,*) "z before reduce"
 534 c        do i=1,dimen
 535 c          write (2,*) i,temp(i)
 536 c        enddo
 537         time00=MPI_Wtime()
 538         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 539      &      MPI_SUM,king,FG_COMM,IERR)
 540         time_reduce=time_reduce+MPI_Wtime()-time00
 541 c        write (2,*) "After REDUCE"
 542 c        call flush(2)
 543       else
 544 #endif
 545 #ifdef TIMING
 546         time01=MPI_Wtime()
 547 #endif
 548         do k=0,2
 549           do i=1,dimen
 550             ind=(i-1)*3+k+1
 551             d_a_tmp(ind)=0.0d0
 552             do j=1,dimen
 553 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 554 c              call flush(2)
 555 c     &         Ginv(i,j),z((j-1)*3+k+1),
 556 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 557               d_a_tmp(ind)=d_a_tmp(ind)
 558      &                         +Ginv(j,i)*z((j-1)*3+k+1)
 559 c              d_a_tmp(ind)=d_a_tmp(ind)
 560 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 561             enddo
 562           enddo
 563         enddo
 564 #ifdef TIMING
 565         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 566 #endif
 567 #ifdef MPI
 568       endif
 569 #endif
 570       return
 571       end
 572 c---------------------------------------------------------------------------
 573 #ifdef GINV_MULT
 574       SUBROUTINE ginv_mult_test(z,d_a_tmp)
 575       include 'DIMENSIONS'
 576       integer dimen
 577 c      include 'COMMON.MD'
 578       double precision z(dimen),d_a_tmp(dimen)
 579       double precision ztmp(dimen/3),dtmp(dimen/3)
 580
 581 c      do i=1,dimen
 582 c        d_a_tmp(i)=0.0d0
 583 c        do j=1,dimen
 584 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 585 c        enddo
 586 c      enddo
 587 c
 588 c      return
 589
 590 !ibm* unroll(3)
 591       do k=0,2
 592        do j=1,dimen/3
 593         ztmp(j)=z((j-1)*3+k+1)
 594        enddo
 595
 596        call alignx(16,ztmp(1))
 597        call alignx(16,dtmp(1))
 598        call alignx(16,Ginv(1,1))
 599
 600        do i=1,dimen/3
 601         dtmp(i)=0.0d0
 602         do j=1,dimen/3
 603            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
 604         enddo
 605        enddo
 606        do i=1,dimen/3
 607         ind=(i-1)*3+k+1
 608         d_a_tmp(ind)=dtmp(i)
 609        enddo
 610       enddo
 611       return
 612       end
 613 #endif
 614 c---------------------------------------------------------------------------
 615       SUBROUTINE fricmat_mult(z,d_a_tmp)
 616       include 'DIMENSIONS'
 617 #ifdef MPI
 618       include 'mpif.h'
 619       integer IERROR
 620 #endif
 621       include 'COMMON.MD'
 622       include 'COMMON.IOUNITS'
 623       include 'COMMON.SETUP'
 624       include 'COMMON.TIME1'
 625 #ifndef LANG0
 626       include 'COMMON.LANGEVIN'
 627 #else
 628       include 'COMMON.LANGEVIN.lang0'
 629 #endif
 630       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 631      &time01
 632 #ifdef MPI
 633       if (nfgtasks.gt.1) then
 634         if (fg_rank.eq.0) then
 635 c The matching BROADCAST for fg processors is called in ERGASTULUM
 636           time00=MPI_Wtime()
 637           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
 638           time_Bcast=time_Bcast+MPI_Wtime()-time00
 639 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
 640         endif
 641 c        call MPI_Barrier(FG_COMM,IERROR)
 642         time00=MPI_Wtime()
 643         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 644      &    MPI_DOUBLE_PRECISION,
 645      &    z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 646 c        write (2,*) "My chunk of z"
 647 c        do i=1,3*my_ng_count
 648 c          write (2,*) i,z(i)
 649 c        enddo
 650         time_scatter=time_scatter+MPI_Wtime()-time00
 651 #ifdef TIMING
 652         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
 653 #endif
 654         time01=MPI_Wtime()
 655         do k=0,2
 656           do i=1,dimen
 657             ind=(i-1)*3+k+1
 658             temp(ind)=0.0d0
 659             do j=1,my_ng_count
 660               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
 661             enddo
 662           enddo
 663         enddo
 664         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 665 c        write (2,*) "Before REDUCE"
 666 c        write (2,*) "d_a_tmp before reduce"
 667 c        do i=1,dimen3
 668 c          write (2,*) i,temp(i)
 669 c        enddo
 670 c        call flush(2)
 671         time00=MPI_Wtime()
 672         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 673      &      MPI_SUM,king,FG_COMM,IERR)
 674         time_reduce=time_reduce+MPI_Wtime()-time00
 675 c        write (2,*) "After REDUCE"
 676 c        call flush(2)
 677       else
 678 #endif
 679 #ifdef TIMING
 680         time01=MPI_Wtime()
 681 #endif
 682         do k=0,2
 683          do i=1,dimen
 684           ind=(i-1)*3+k+1
 685           d_a_tmp(ind)=0.0d0
 686           do j=1,dimen
 687              d_a_tmp(ind)=d_a_tmp(ind)
 688      &                           -fricmat(j,i)*z((j-1)*3+k+1)
 689           enddo
 690          enddo
 691         enddo
 692 #ifdef TIMING
 693         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 694 #endif
 695 #ifdef MPI
 696       endif
 697 #endif
 698 c      write (iout,*) "Vector d_a"
 699 c      do i=1,dimen3
 700 c        write (2,*) i,d_a_tmp(i)
 701 c      enddo
 702       return
 703       end