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