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