source/unres/src_MD/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,zcopy(dimen3)
 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      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 524          do i=1,3*my_ng_count
 525            z(i)=zcopy(i)
 526          enddo
 527 c        write (2,*) "My chunk of z"
 528 c        do i=1,3*my_ng_count
 529 c          write (2,*) i,z(i)
 530 c        enddo
 531 c        write (2,*) "After SCATTERV"
 532 c        call flush(2)
 533 c        write (2,*) "MPI_Wtime",MPI_Wtime()
 534         time_scatter=time_scatter+MPI_Wtime()-time00
 535 #ifdef TIMING
 536         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 537 #endif
 538 c        write (2,*) "time_scatter",time_scatter
 539 c        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 540 c     &    my_ng_count
 541 c        call flush(2)
 542         time01=MPI_Wtime()
 543         do k=0,2
 544           do i=1,dimen
 545             ind=(i-1)*3+k+1
 546             temp(ind)=0.0d0
 547             do j=1,my_ng_count
 548 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 549 c     &         Ginv(i,j),z((j-1)*3+k+1),
 550 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 551 c              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 552               temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1)
 553             enddo
 554           enddo
 555         enddo
 556         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 557 c        write (2,*) "Before REDUCE"
 558 c        call flush(2)
 559 c        write (2,*) "z before reduce"
 560 c        do i=1,dimen
 561 c          write (2,*) i,temp(i)
 562 c        enddo
 563         time00=MPI_Wtime()
 564         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 565      &      MPI_SUM,king,FG_COMM,IERR)
 566         time_reduce=time_reduce+MPI_Wtime()-time00
 567 c        write (2,*) "After REDUCE"
 568 c        call flush(2)
 569       else
 570 #endif
 571 #ifdef TIMING
 572         time01=MPI_Wtime()
 573 #endif
 574         do k=0,2
 575           do i=1,dimen
 576             ind=(i-1)*3+k+1
 577             d_a_tmp(ind)=0.0d0
 578             do j=1,dimen
 579 c              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 580 c              call flush(2)
 581 c     &         Ginv(i,j),z((j-1)*3+k+1),
 582 c     &          Ginv(i,j)*z((j-1)*3+k+1)
 583               d_a_tmp(ind)=d_a_tmp(ind)
 584      &                         +Ginv(j,i)*z((j-1)*3+k+1)
 585 c              d_a_tmp(ind)=d_a_tmp(ind)
 586 c     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 587             enddo
 588           enddo
 589         enddo
 590 #ifdef TIMING
 591         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 592 #endif
 593 #ifdef MPI
 594       endif
 595 #endif
 596       return
 597       end
 598 c---------------------------------------------------------------------------
 599 #ifdef GINV_MULT
 600       SUBROUTINE ginv_mult_test(z,d_a_tmp)
 601       include 'DIMENSIONS'
 602       integer dimen
 603 c      include 'COMMON.MD'
 604       double precision z(dimen),d_a_tmp(dimen)
 605       double precision ztmp(dimen/3),dtmp(dimen/3)
 606
 607 c      do i=1,dimen
 608 c        d_a_tmp(i)=0.0d0
 609 c        do j=1,dimen
 610 c          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 611 c        enddo
 612 c      enddo
 613 c
 614 c      return
 615
 616 !ibm* unroll(3)
 617       do k=0,2
 618        do j=1,dimen/3
 619         ztmp(j)=z((j-1)*3+k+1)
 620        enddo
 621
 622        call alignx(16,ztmp(1))
 623        call alignx(16,dtmp(1))
 624        call alignx(16,Ginv(1,1))
 625
 626        do i=1,dimen/3
 627         dtmp(i)=0.0d0
 628         do j=1,dimen/3
 629            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
 630         enddo
 631        enddo
 632        do i=1,dimen/3
 633         ind=(i-1)*3+k+1
 634         d_a_tmp(ind)=dtmp(i)
 635        enddo
 636       enddo
 637       return
 638       end
 639 #endif
 640 c---------------------------------------------------------------------------
 641       SUBROUTINE fricmat_mult(z,d_a_tmp)
 642       include 'DIMENSIONS'
 643 #ifdef MPI
 644       include 'mpif.h'
 645       integer IERROR
 646 #endif
 647       include 'COMMON.MD'
 648       include 'COMMON.IOUNITS'
 649       include 'COMMON.SETUP'
 650       include 'COMMON.TIME1'
 651 #ifndef LANG0
 652       include 'COMMON.LANGEVIN'
 653 #else
 654       include 'COMMON.LANGEVIN.lang0'
 655 #endif
 656       double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00
 657      &,time01,zcopy(dimen3)
 658 #ifdef MPI
 659       if (nfgtasks.gt.1) then
 660         if (fg_rank.eq.0) then
 661 c The matching BROADCAST for fg processors is called in ERGASTULUM
 662           time00=MPI_Wtime()
 663           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
 664           time_Bcast=time_Bcast+MPI_Wtime()-time00
 665 c          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
 666         endif
 667 c        call MPI_Barrier(FG_COMM,IERROR)
 668         time00=MPI_Wtime()
 669         call MPI_Scatterv(z,ng_counts(0),ng_start(0),
 670      &    MPI_DOUBLE_PRECISION,
 671      &    zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 672
 673          do i=1,3*my_ng_count
 674           z(i)=zcopy(i)
 675          enddo
 676 c        write (2,*) "My chunk of z"
 677 c        do i=1,3*my_ng_count
 678 c          write (2,*) i,z(i)
 679 c        enddo
 680         time_scatter=time_scatter+MPI_Wtime()-time00
 681 #ifdef TIMING
 682         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
 683 #endif
 684         time01=MPI_Wtime()
 685         do k=0,2
 686           do i=1,dimen
 687             ind=(i-1)*3+k+1
 688             temp(ind)=0.0d0
 689             do j=1,my_ng_count
 690               temp(ind)=temp(ind)-fricmat(j,i)*z((j-1)*3+k+1)
 691             enddo
 692           enddo
 693         enddo
 694         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 695 c        write (2,*) "Before REDUCE"
 696 c        write (2,*) "d_a_tmp before reduce"
 697 c        do i=1,dimen3
 698 c          write (2,*) i,temp(i)
 699 c        enddo
 700 c        call flush(2)
 701         time00=MPI_Wtime()
 702         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,
 703      &      MPI_SUM,king,FG_COMM,IERR)
 704         time_reduce=time_reduce+MPI_Wtime()-time00
 705 c        write (2,*) "After REDUCE"
 706 c        call flush(2)
 707       else
 708 #endif
 709 #ifdef TIMING
 710         time01=MPI_Wtime()
 711 #endif
 712         do k=0,2
 713          do i=1,dimen
 714           ind=(i-1)*3+k+1
 715           d_a_tmp(ind)=0.0d0
 716           do j=1,dimen
 717              d_a_tmp(ind)=d_a_tmp(ind)
 718      &                           -fricmat(j,i)*z((j-1)*3+k+1)
 719           enddo
 720          enddo
 721         enddo
 722 #ifdef TIMING
 723         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 724 #endif
 725 #ifdef MPI
 726       endif
 727 #endif
 728 c      write (iout,*) "Vector d_a"
 729 c      do i=1,dimen3
 730 c        write (2,*) i,d_a_tmp(i)
 731 c      enddo
 732       return
 733       end