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        write (iout,*) "The number of degrees of freedom ",dimen3
 186 c  Zeroing out A and fricmat
 187       do i=1,dimen
 188         do j=1,dimen
 189           A(i,j)=0.0D0
 190         enddo
 191       enddo
 192 c  Diagonal elements of the dC part of A and the respective friction coefficients
 193       ind=1
 194       ind1=0
 195       do i=nnt,nct-1
 196         ind=ind+1
 197         ind1=ind1+1
 198         coeff=0.25d0*IP
 199         massvec(ind1)=mp
 200         Gmat(ind,ind)=coeff
 201         A(ind1,ind)=0.5d0
 202       enddo
 203
 204 c  Off-diagonal elements of the dC part of A
 205       k=3
 206       do i=1,nct-nnt
 207         do j=1,i
 208           A(i,j)=1.0d0
 209         enddo
 210       enddo
 211 c  Diagonal elements of the dX part of A and the respective friction coefficients
 212       m=nct-nnt
 213       m1=nct-nnt+1
 214       ind=0
 215       ind1=0
 216       msc(21)=1.0d0
 217       do i=nnt,nct
 218         ind=ind+1
 219         ii = ind+m
 220         iti=itype(i)
 221         massvec(ii)=msc(iabs(iti))
 222         if (iti.ne.10 .and. iti.ne.ntyp1) then
 223           ind1=ind1+1
 224           ii1= ind1+m1
 225           A(ii,ii1)=1.0d0
 226           Gmat(ii1,ii1)=ISC(iabs(iti))
 227         endif
 228       enddo
 229 c  Off-diagonal elements of the dX part of A
 230       ind=0
 231       k=nct-nnt
 232       do i=nnt,nct
 233         iti=itype(i)
 234         ind=ind+1
 235         do j=nnt,i
 236           ii = ind
 237           jj = j-nnt+1
 238           A(k+ii,jj)=1.0d0
 239         enddo
 240       enddo
 241       if (lprn) then
 242         write (iout,*)
 243         write (iout,*) "Vector massvec"
 244         do i=1,dimen1
 245           write (iout,*) i,massvec(i)
 246         enddo
 247         write (iout,'(//a)') "A"
 248         call matout(dimen,dimen1,maxres2,maxres2,A)
 249       endif
 250
 251 c Calculate the G matrix (store in Gmat)
 252       do k=1,dimen
 253        do i=1,dimen
 254          dtdi=0.0d0
 255          do j=1,dimen1
 256            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 257          enddo
 258          Gmat(k,i)=Gmat(k,i)+dtdi
 259        enddo
 260       enddo
 261
 262       if (lprn) then
 263         write (iout,'(//a)') "Gmat"
 264         call matout(dimen,dimen,maxres2,maxres2,Gmat)
 265       endif
 266       do i=1,dimen
 267         do j=1,dimen
 268           Ginv(i,j)=0.0d0
 269           Gcopy(i,j)=Gmat(i,j)
 270         enddo
 271         Ginv(i,i)=1.0d0
 272       enddo
 273 c Invert the G matrix
 274       call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
 275       if (lprn) then
 276         write (iout,'(//a)') "Ginv"
 277         call matout(dimen,dimen,maxres2,maxres2,Ginv)
 278       endif
 279 #ifdef MPI
 280       if (nfgtasks.gt.1) then
 281         myginv_ng_count=maxres2*my_ng_count
 282         call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
 283      &    nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 284         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
 285      &    nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 286         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 287         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 288         call flush(iout)
 289 c        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 290 c     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 291 c     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 292 c        call MPI_Barrier(FG_COMM,IERR)
 293         time00=MPI_Wtime()
 294         call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 295      &    nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
 296      &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 297 #ifdef TIMING
 298         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 299 #endif
 300         do i=1,dimen
 301           do j=1,2*my_ng_count
 302             ginv(j,i)=gcopy(i,j)
 303           enddo
 304         enddo
 305 c        write (iout,*) "Master's chunk of ginv"
 306 c        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 307       endif
 308 #endif
 309       if (osob) then
 310         write (iout,*) "The G matrix is singular."
 311         stop
 312       endif
 313 c Compute G**(-1/2) and G**(1/2)
 314       ind=0
 315       do i=1,dimen
 316         do j=1,i
 317           ind=ind+1
 318           Ghalf(ind)=Gmat(i,j)
 319         enddo
 320       enddo
 321       call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
 322      &  ierr,iwork)
 323       if (lprn) then
 324         write (iout,'(//a)')
 325      &   "Eigenvectors and eigenvalues of the G matrix"
 326         call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
 327       endif
 328       do i=1,dimen
 329         sqreig(i)=dsqrt(Geigen(i))
 330       enddo
 331       do i=1,dimen
 332         do j=1,dimen
 333           Gsqrp(i,j)=0.0d0
 334           Gsqrm(i,j)=0.0d0
 335           Gcopy(i,j)=0.0d0
 336           do k=1,dimen
 337             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 338             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 339             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 340           enddo
 341         enddo
 342       enddo
 343       if (lprn) then
 344         write (iout,*) "Comparison of original and restored G"
 345         do i=1,dimen
 346           do j=1,dimen
 347             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),
 348      &        Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 349           enddo
 350         enddo
 351       endif
 352       return
 353       end
 354 c-------------------------------------------------------------------------------
 355       SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
 356       implicit real*8 (a-h,o-z)
 357       include 'DIMENSIONS'
 358       include 'COMMON.IOUNITS'
 359       double precision A(LM2,LM3),B(LM2)
 360       KA=1
 361       KC=6
 362     1 KB=MIN0(KC,NC)
 363       WRITE(IOUT,600) (I,I=KA,KB)
 364       WRITE(IOUT,601) (B(I),I=KA,KB)
 365       WRITE(IOUT,602)
 366     2 N=0
 367       DO 3  I=1,NR
 368       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 369       N=N+1
 370       IF(N.LT.10) GO TO 3
 371       WRITE(IOUT,602)
 372       N=0
 373     3 CONTINUE
 374     4 IF (KB.EQ.NC) RETURN
 375       KA=KC+1
 376       KC=KC+6
 377       GO TO 1
 378   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 379   601 FORMAT (/5X,10(1PE11.4))
 380   602 FORMAT (2H  )
 381   603 FORMAT (I5,10F11.5)
 382   604 FORMAT (1H1)
 383       END
 384 c-------------------------------------------------------------------------------
 385       SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
 386       implicit real*8 (a-h,o-z)
 387       include 'DIMENSIONS'
 388       include 'COMMON.IOUNITS'
 389       double precision A(LM2,LM3)
 390       KA=1
 391       KC=6
 392     1 KB=MIN0(KC,NC)
 393       WRITE(IOUT,600) (I,I=KA,KB)
 394       WRITE(IOUT,602)
 395     2 N=0
 396       DO 3  I=1,NR
 397       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 398       N=N+1
 399       IF(N.LT.10) GO TO 3
 400       WRITE(IOUT,602)
 401       N=0
 402     3 CONTINUE
 403     4 IF (KB.EQ.NC) RETURN
 404       KA=KC+1
 405       KC=KC+6
 406       GO TO 1
 407   600 FORMAT (//5x,9I11)
 408   602 FORMAT (2H  )
 409   603 FORMAT (I5,10F11.3)
 410   604 FORMAT (1H1)
 411       END
 412 c-------------------------------------------------------------------------------
 413       SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
 414       implicit real*8 (a-h,o-z)
 415       include 'DIMENSIONS'
 416       include 'COMMON.IOUNITS'
 417       double precision A(LM2,LM3)
 418       KA=1
 419       KC=21
 420     1 KB=MIN0(KC,NC)
 421       WRITE(IOUT,600) (I,I=KA,KB)
 422       WRITE(IOUT,602)
 423     2 N=0
 424       DO 3  I=1,NR
 425       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 426       N=N+1
 427       IF(N.LT.3) GO TO 3
 428       WRITE(IOUT,602)
 429       N=0
 430     3 CONTINUE
 431     4 IF (KB.EQ.NC) RETURN
 432       KA=KC+1
 433       KC=KC+21
 434       GO TO 1
 435   600 FORMAT (//5x,7(3I5,2x))
 436   602 FORMAT (2H  )
 437   603 FORMAT (I5,7(3F5.1,2x))
 438   604 FORMAT (1H1)
 439       END
 440 c-------------------------------------------------------------------------------
 441       SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
 442       implicit real*8 (a-h,o-z)
 443       include 'DIMENSIONS'
 444       include 'COMMON.IOUNITS'
 445       double precision A(LM2,LM3)
 446       KA=1
 447       KC=12
 448     1 KB=MIN0(KC,NC)
 449       WRITE(IOUT,600) (I,I=KA,KB)
 450       WRITE(IOUT,602)
 451     2 N=0
 452       DO 3  I=1,NR
 453       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 454       N=N+1
 455       IF(N.LT.3) GO TO 3
 456       WRITE(IOUT,602)
 457       N=0
 458     3 CONTINUE
 459     4 IF (KB.EQ.NC) RETURN
 460       KA=KC+1
 461       KC=KC+12
 462       GO TO 1
 463   600 FORMAT (//5x,4(3I9,2x))
 464   602 FORMAT (2H  )
 465   603 FORMAT (I5,4(3F9.3,2x))
 466   604 FORMAT (1H1)
 467       END
 468 c---------------------------------------------------------------------------
 469       SUBROUTINE ginv_mult(z,d_a_tmp)
 470       implicit real*8 (a-h,o-z)
 471       include 'DIMENSIONS'
 472 #ifdef MPI
 473       include 'mpif.h'
 474       integer ierr
 475 #endif
 476       include 'COMMON.SETUP'
 477       include 'COMMON.TIME1'
 478       include 'COMMON.MD'
 479       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 480      &,time01,zcopy(dimen3)
 481 #ifdef MPI
 482       if (nfgtasks.gt.1) then
 483         if (fg_rank.eq.0) then
 484 c The matching BROADCAST for fg processors is called in ERGASTULUM
 485           time00=MPI_Wtime()
 486           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
 487           time_Bcast=time_Bcast+MPI_Wtime()-time00
 488 c          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
 489         endif
 490 c        write (2,*) "time00",time00
 491 c        write (2,*) "Before Scatterv"
 492 c        call flush(2)
 493 c        write (2,*) "Whole z (for FG master)"
 494 c        do i=1,dimen
 495 c          write (2,*) i,z(i)
 496 c        enddo
 497 c        call MPI_Barrier(FG_COMM,IERROR)
 498         time00=MPI_Wtime()
 499         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 500      &    MPI_DOUBLE_PRECISION,
 501      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 502          do i=1,3*my_ng_count
 503            z(i)=zcopy(i)
 504          enddo
 505 c        write (2,*) "My chunk of z"
 506 c        do i=1,3*my_ng_count
 507 c          write (2,*) i,z(i)
 508 c        enddo
 509 c        write (2,*) "After SCATTERV"
 510 c        call flush(2)
 511 c        write (2,*) "MPI_Wtime",MPI_Wtime()
 512         time_scatter=time_scatter+MPI_Wtime()-time00
 513 #ifdef TIMING
 514         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 515 #endif
 516 c        write (2,*) "time_scatter",time_scatter
 517 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 518 c     &    my_ng_count
 519 c        call flush(2)
 520         time01=MPI_Wtime()
 521         do k=0,2
 522           do i=1,dimen
 523             ind=(i-1)*3+k+1
 524             temp(ind)=0.0d0
 525             do j=1,my_ng_count
 526 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 527 c     &         Ginv(i,j),z((j-1)*3+k+1),
 528 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 529 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 530               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
 531             enddo
 532           enddo
 533         enddo
 534         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 535 c        write (2,*) "Before REDUCE"
 536 c        call flush(2)
 537 c        write (2,*) "z before reduce"
 538 c        do i=1,dimen
 539 c          write (2,*) i,temp(i)
 540 c        enddo
 541         time00=MPI_Wtime()
 542         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 543      &      MPI_SUM,king,FG_COMM,IERR)
 544         time_reduce=time_reduce+MPI_Wtime()-time00
 545 c        write (2,*) "After REDUCE"
 546 c        call flush(2)
 547       else
 548 #endif
 549 #ifdef TIMING
 550         time01=MPI_Wtime()
 551 #endif
 552         do k=0,2
 553           do i=1,dimen
 554             ind=(i-1)*3+k+1
 555             d_a_tmp(ind)=0.0d0
 556             do j=1,dimen
 557 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 558 c              call flush(2)
 559 c     &         Ginv(i,j),z((j-1)*3+k+1),
 560 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 561               d_a_tmp(ind)=d_a_tmp(ind)
 562      &                         +Ginv(j,i)*z((j-1)*3+k+1)
 563 c              d_a_tmp(ind)=d_a_tmp(ind)
 564 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 565             enddo
 566           enddo
 567         enddo
 568 #ifdef TIMING
 569         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 570 #endif
 571 #ifdef MPI
 572       endif
 573 #endif
 574       return
 575       end
 576 c---------------------------------------------------------------------------
 577 #ifdef GINV_MULT
 578       SUBROUTINE ginv_mult_test(z,d_a_tmp)
 579       include 'DIMENSIONS'
 580       integer dimen
 581 c      include 'COMMON.MD'
 582       double precision z(dimen),d_a_tmp(dimen)
 583       double precision ztmp(dimen/3),dtmp(dimen/3)
 584
 585 c      do i=1,dimen
 586 c        d_a_tmp(i)=0.0d0
 587 c        do j=1,dimen
 588 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 589 c        enddo
 590 c      enddo
 591 c
 592 c      return
 593
 594 !ibm* unroll(3)
 595       do k=0,2
 596        do j=1,dimen/3
 597         ztmp(j)=z((j-1)*3+k+1)
 598        enddo
 599
 600        call alignx(16,ztmp(1))
 601        call alignx(16,dtmp(1))
 602        call alignx(16,Ginv(1,1))
 603
 604        do i=1,dimen/3
 605         dtmp(i)=0.0d0
 606         do j=1,dimen/3
 607            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
 608         enddo
 609        enddo
 610        do i=1,dimen/3
 611         ind=(i-1)*3+k+1
 612         d_a_tmp(ind)=dtmp(i)
 613        enddo
 614       enddo
 615       return
 616       end
 617 #endif
 618 c---------------------------------------------------------------------------
 619       SUBROUTINE fricmat_mult(z,d_a_tmp)
 620       include 'DIMENSIONS'
 621 #ifdef MPI
 622       include 'mpif.h'
 623       integer IERROR
 624 #endif
 625       include 'COMMON.MD'
 626       include 'COMMON.IOUNITS'
 627       include 'COMMON.SETUP'
 628       include 'COMMON.TIME1'
 629 #ifndef LANG0
 630       include 'COMMON.LANGEVIN'
 631 #else
 632       include 'COMMON.LANGEVIN.lang0'
 633 #endif
 634       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 635      &,time01,zcopy(dimen3)
 636 #ifdef MPI
 637       if (nfgtasks.gt.1) then
 638         if (fg_rank.eq.0) then
 639 c The matching BROADCAST for fg processors is called in ERGASTULUM
 640           time00=MPI_Wtime()
 641           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
 642           time_Bcast=time_Bcast+MPI_Wtime()-time00
 643 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
 644         endif
 645 c        call MPI_Barrier(FG_COMM,IERROR)
 646         time00=MPI_Wtime()
 647         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 648      &    MPI_DOUBLE_PRECISION,
 649      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 650
 651          do i=1,3*my_ng_count
 652           z(i)=zcopy(i)
 653          enddo
 654 c        write (2,*) "My chunk of z"
 655 c        do i=1,3*my_ng_count
 656 c          write (2,*) i,z(i)
 657 c        enddo
 658         time_scatter=time_scatter+MPI_Wtime()-time00
 659 #ifdef TIMING
 660         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
 661 #endif
 662         time01=MPI_Wtime()
 663         do k=0,2
 664           do i=1,dimen
 665             ind=(i-1)*3+k+1
 666             temp(ind)=0.0d0
 667             do j=1,my_ng_count
 668               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
 669             enddo
 670           enddo
 671         enddo
 672         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 673 c        write (2,*) "Before REDUCE"
 674 c        write (2,*) "d_a_tmp before reduce"
 675 c        do i=1,dimen3
 676 c          write (2,*) i,temp(i)
 677 c        enddo
 678 c        call flush(2)
 679         time00=MPI_Wtime()
 680         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 681      &      MPI_SUM,king,FG_COMM,IERR)
 682         time_reduce=time_reduce+MPI_Wtime()-time00
 683 c        write (2,*) "After REDUCE"
 684 c        call flush(2)
 685       else
 686 #endif
 687 #ifdef TIMING
 688         time01=MPI_Wtime()
 689 #endif
 690         do k=0,2
 691          do i=1,dimen
 692           ind=(i-1)*3+k+1
 693           d_a_tmp(ind)=0.0d0
 694           do j=1,dimen
 695              d_a_tmp(ind)=d_a_tmp(ind)
 696      &                           -fricmat(j,i)*z((j-1)*3+k+1)
 697           enddo
 698          enddo
 699         enddo
 700 #ifdef TIMING
 701         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 702 #endif
 703 #ifdef MPI
 704       endif
 705 #endif
 706 c      write (iout,*) "Vector d_a"
 707 c      do i=1,dimen3
 708 c        write (2,*) i,d_a_tmp(i)
 709 c      enddo
 710       return
 711       end