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