source/unres/REMD.f90

   1       module REMD
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use MD_data
   5       use REMD_data
   6       use muca_md
   7
   8       implicit none
   9 !-----------------------------------------------------------------------------
  10 !
  11 !
  12 !-----------------------------------------------------------------------------
  13       contains
  14 !-----------------------------------------------------------------------------
  15 ! lagrangian_lesyng.F
  16 !-----------------------------------------------------------------------------
  17       subroutine lagrangian
  18 !-------------------------------------------------------------------------
  19 !  This subroutine contains the total lagrangain from which the accelerations
  20 !  are obtained.  For numerical gradient checking, the derivetive of the
  21 !  lagrangian in the velocities and coordinates are calculated seperately
  22 !-------------------------------------------------------------------------
  23 !       implicit real*8 (a-h,o-z)
  24 !       include 'DIMENSIONS'
  25       use comm_cipiszcze
  26       use energy_data
  27       use geometry_data, only: nres
  28       use control_data
  29 #ifdef MPI
  30        include 'mpif.h'
  31       real(kind=8) :: time00
  32 #endif
  33 !       include 'COMMON.VAR'
  34 !       include 'COMMON.CHAIN'
  35 !       include 'COMMON.DERIV'
  36 !       include 'COMMON.GEO'
  37 !       include 'COMMON.LOCAL'
  38 !       include 'COMMON.INTERACT'
  39 !       include 'COMMON.MD'
  40 !       include 'COMMON.IOUNITS'
  41 !       include 'COMMON.CONTROL'
  42 !       include 'COMMON.MUCA'
  43 !      include 'COMMON.TIME1'
  44
  45        integer :: i,j,ind,itime
  46        real(kind=8) :: zapas(6*nres) !,muca_factor      !maxres6=6*maxres
  47        logical :: lprn = .false.
  48 !el       common /cipiszcze/ itime
  49        itime = itt_comm
  50
  51 #ifdef TIMING
  52        time00=MPI_Wtime()
  53 #endif
  54        do j=1,3
  55          zapas(j)=-gcart(j,0)
  56        enddo
  57       ind=3
  58       if (lprn) then
  59         write (iout,*) "Potential forces backbone"
  60       endif
  61       do i=nnt,nct-1
  62         if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
  63           i,(-gcart(j,i),j=1,3)
  64         do j=1,3
  65           ind=ind+1
  66           zapas(ind)=-gcart(j,i)
  67         enddo
  68       enddo
  69       if (lprn) write (iout,*) "Potential forces sidechain"
  70       do i=nnt,nct
  71         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
  72           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
  73              i,(-gcart(j,i),j=1,3)
  74           do j=1,3
  75             ind=ind+1
  76             zapas(ind)=-gxcart(j,i)
  77           enddo
  78         endif
  79       enddo
  80
  81       call ginv_mult(zapas,d_a_work)
  82
  83       do j=1,3
  84         d_a(j,0)=d_a_work(j)
  85       enddo
  86       ind=3
  87       do i=nnt,nct-1
  88         do j=1,3
  89           ind=ind+1
  90           d_a(j,i)=d_a_work(ind)
  91         enddo
  92       enddo
  93       do i=nnt,nct
  94         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
  95           do j=1,3
  96             ind=ind+1
  97             d_a(j,i+nres)=d_a_work(ind)
  98           enddo
  99         endif
 100       enddo
 101
 102       if(lmuca) then
 103        imtime=imtime+1
 104        if(mucadyn.gt.0) call muca_update(potE)
 105        factor=muca_factor(potE)*t_bath*Rb
 106
 107 !d       print *,'lmuca ',factor,potE
 108        do j=1,3
 109           d_a(j,0)=d_a(j,0)*factor
 110        enddo
 111        do i=nnt,nct-1
 112          do j=1,3
 113           d_a(j,i)=d_a(j,i)*factor
 114          enddo
 115        enddo
 116        do i=nnt,nct
 117          do j=1,3
 118           d_a(j,i+nres)=d_a(j,i+nres)*factor
 119          enddo
 120        enddo
 121
 122       endif
 123
 124       if (lprn) then
 125         write(iout,*) 'acceleration 3D'
 126         write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
 127         do i=nnt,nct-1
 128          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
 129         enddo
 130         do i=nnt,nct
 131          write (iout,'(i3,3f10.5,3x,3f10.5)') &
 132            i+nres,(d_a(j,i+nres),j=1,3)
 133         enddo
 134       endif
 135 #ifdef TIMING
 136       time_lagrangian=time_lagrangian+MPI_Wtime()-time00
 137 #endif
 138       return
 139       end subroutine lagrangian
 140 !-----------------------------------------------------------------------------
 141       subroutine setup_MD_matrices
 142
 143       use geometry_data, only: nres,nside
 144       use control_data
 145       use MPI_data
 146       use energy_data
 147       use geometry, only:int_bounds
 148       use md_calc
 149 !      implicit real*8 (a-h,o-z)
 150 !      include 'DIMENSIONS'
 151 #ifdef MPI
 152       include 'mpif.h'
 153       integer :: ierror
 154       real(kind=8) :: time00
 155 #endif
 156 !      include 'COMMON.SETUP'
 157 !      include 'COMMON.VAR'
 158 !      include 'COMMON.CHAIN'
 159 !      include 'COMMON.DERIV'
 160 !      include 'COMMON.GEO'
 161 !      include 'COMMON.LOCAL'
 162 !      include 'COMMON.INTERACT'
 163 !      include 'COMMON.MD'
 164 !#ifndef LANG0
 165 !      include 'COMMON.LANGEVIN'
 166 !#else
 167 !      include 'COMMON.LANGEVIN.lang0'
 168 !#endif
 169 !      include 'COMMON.IOUNITS'
 170 !      include 'COMMON.TIME1'
 171       logical :: lprn = .false.
 172       logical :: osob
 173       real(kind=8) :: dtdi
 174       real(kind=8),dimension(2*nres) :: massvec,sqreig  !(maxres2) maxres2=2*maxres
 175 !el      real(kind=8),dimension(:),allocatable :: Ghalf
 176 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf   !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
 177 !el      real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2)
 178 !el      real(kind=8),dimension(:,:),allocatable :: Gcopy
 179       real(kind=8),dimension(8*6*nres) :: work  !(8*maxres6)
 180       integer,dimension(6*nres) :: iwork        !(maxres6) maxres6=6*maxres
 181 !el      common /przechowalnia/ Gcopy,Ghalf
 182       real(kind=8) :: coeff
 183       integer :: i,j,ind,ind1,k,ii,jj,m,m1,ii1,iti,nres2,ierr
 184       nres2=2*nres
 185
 186       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
 187       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
 188 !
 189 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 190 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
 191 !
 192 ! Determine the number of degrees of freedom (dimen) and the number of
 193 ! sites (dimen1)
 194       dimen=(nct-nnt+1)+nside
 195       dimen1=(nct-nnt)+(nct-nnt+1)
 196       dimen3=dimen*3
 197 #ifdef MPI
 198       if (nfgtasks.gt.1) then
 199       time00=MPI_Wtime()
 200       call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR)
 201       time_Bcast=time_Bcast+MPI_Wtime()-time00
 202       call int_bounds(dimen,igmult_start,igmult_end)
 203       igmult_start=igmult_start-1
 204       call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,&
 205           ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 206       my_ng_count=igmult_end-igmult_start
 207       call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,&
 208           MPI_INTEGER,FG_COMM,IERROR)
 209       write (iout,*) 'Processor:',fg_rank,' CG group',kolor,&
 210        ' absolute rank',myrank,' igmult_start',igmult_start,&
 211        ' igmult_end',igmult_end,' count',my_ng_count
 212       write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
 213       write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
 214       call flush(iout)
 215       else
 216 #endif
 217       igmult_start=1
 218       igmult_end=dimen
 219       my_ng_count=dimen
 220 #ifdef MPI
 221       endif
 222 #endif
 223 !      write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3
 224 !  Zeroing out A and fricmat
 225       do i=1,dimen
 226         do j=1,dimen
 227           A(i,j)=0.0D0
 228         enddo
 229       enddo
 230 !  Diagonal elements of the dC part of A and the respective friction coefficients
 231       ind=1
 232       ind1=0
 233       do i=nnt,nct-1
 234         ind=ind+1
 235         ind1=ind1+1
 236         coeff=0.25d0*IP
 237         massvec(ind1)=mp
 238         Gmat(ind,ind)=coeff
 239         A(ind1,ind)=0.5d0
 240       enddo
 241
 242 !  Off-diagonal elements of the dC part of A
 243       k=3
 244       do i=1,nct-nnt
 245         do j=1,i
 246           A(i,j)=1.0d0
 247         enddo
 248       enddo
 249 !  Diagonal elements of the dX part of A and the respective friction coefficients
 250       m=nct-nnt
 251       m1=nct-nnt+1
 252       ind=0
 253       ind1=0
 254       msc(ntyp1)=1.0d0
 255       do i=nnt,nct
 256         ind=ind+1
 257         ii = ind+m
 258         iti=itype(i)
 259         massvec(ii)=msc(iabs(iti))
 260         if (iti.ne.10 .and. iti.ne.ntyp1) then
 261           ind1=ind1+1
 262           ii1= ind1+m1
 263           A(ii,ii1)=1.0d0
 264           Gmat(ii1,ii1)=ISC(iabs(iti))
 265         endif
 266       enddo
 267 !  Off-diagonal elements of the dX part of A
 268       ind=0
 269       k=nct-nnt
 270       do i=nnt,nct
 271         iti=itype(i)
 272         ind=ind+1
 273         do j=nnt,i
 274           ii = ind
 275           jj = j-nnt+1
 276           A(k+ii,jj)=1.0d0
 277         enddo
 278       enddo
 279       if (lprn) then
 280         write (iout,*)
 281         write (iout,*) "Vector massvec"
 282         do i=1,dimen1
 283           write (iout,*) i,massvec(i)
 284         enddo
 285         write (iout,'(//a)') "A"
 286         call matout(dimen,dimen1,nres2,nres2,A)
 287       endif
 288
 289 ! Calculate the G matrix (store in Gmat)
 290       do k=1,dimen
 291        do i=1,dimen
 292          dtdi=0.0d0
 293          do j=1,dimen1
 294            dtdi=dtdi+A(j,k)*A(j,i)*massvec(j)
 295          enddo
 296          Gmat(k,i)=Gmat(k,i)+dtdi
 297        enddo
 298       enddo
 299
 300       if (lprn) then
 301         write (iout,'(//a)') "Gmat"
 302         call matout(dimen,dimen,nres2,nres2,Gmat)
 303       endif
 304       do i=1,dimen
 305         do j=1,dimen
 306           Ginv(i,j)=0.0d0
 307           Gcopy(i,j)=Gmat(i,j)
 308         enddo
 309         Ginv(i,i)=1.0d0
 310       enddo
 311 ! Invert the G matrix
 312       call MATINVERT(dimen,nres2,Gcopy,Ginv,osob)
 313       if (lprn) then
 314         write (iout,'(//a)') "Ginv"
 315         call matout(dimen,dimen,nres2,nres2,Ginv)
 316       endif
 317 #ifdef MPI
 318       if (nfgtasks.gt.1) then
 319         myginv_ng_count=nres2*my_ng_count
 320         call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,&
 321           nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
 322         call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,&
 323           nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
 324         write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
 325         write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
 326         call flush(iout)
 327 !        call MPI_Scatterv(ginv(1,1),nginv_counts(0),
 328 !     &    nginv_start(0),MPI_DOUBLE_PRECISION,ginv,
 329 !     &    myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 330 !        call MPI_Barrier(FG_COMM,IERR)
 331         time00=MPI_Wtime()
 332         call MPI_Scatterv(ginv(1,1),nginv_counts(0),&
 333           nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),&
 334           myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 335 #ifdef TIMING
 336         time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
 337 #endif
 338         do i=1,dimen
 339           do j=1,2*my_ng_count
 340             ginv(j,i)=gcopy(i,j)
 341           enddo
 342         enddo
 343 !        write (iout,*) "Master's chunk of ginv"
 344 !        call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv)
 345       endif
 346 #endif
 347       if (osob) then
 348         write (iout,*) "The G matrix is singular."
 349         stop
 350       endif
 351 ! Compute G**(-1/2) and G**(1/2)
 352       ind=0
 353       do i=1,dimen
 354         do j=1,i
 355           ind=ind+1
 356           Ghalf(ind)=Gmat(i,j)
 357         enddo
 358       enddo
 359       call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,&
 360         ierr,iwork)
 361       if (lprn) then
 362         write (iout,'(//a)') &
 363          "Eigenvectors and eigenvalues of the G matrix"
 364         call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen)
 365       endif
 366       do i=1,dimen
 367         sqreig(i)=dsqrt(Geigen(i))
 368       enddo
 369       do i=1,dimen
 370         do j=1,dimen
 371           Gsqrp(i,j)=0.0d0
 372           Gsqrm(i,j)=0.0d0
 373           Gcopy(i,j)=0.0d0
 374           do k=1,dimen
 375             Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k)
 376             Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k)
 377             Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k)
 378           enddo
 379         enddo
 380       enddo
 381       if (lprn) then
 382         write (iout,*) "Comparison of original and restored G"
 383         do i=1,dimen
 384           do j=1,dimen
 385             write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),&
 386               Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j)
 387           enddo
 388         enddo
 389       endif
 390       deallocate(Gcopy)
 391       return
 392       end subroutine setup_MD_matrices
 393 !-----------------------------------------------------------------------------
 394       subroutine EIGOUT(NC,NR,LM2,LM3,A,B)
 395
 396 !      implicit real*8 (a-h,o-z)
 397 !      include 'DIMENSIONS'
 398 !      include 'COMMON.IOUNITS'
 399       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 400       real(kind=8) :: A(LM2,LM3),B(LM2)
 401       KA=1
 402       KC=6
 403     1 KB=MIN0(KC,NC)
 404       WRITE(IOUT,600) (I,I=KA,KB)
 405       WRITE(IOUT,601) (B(I),I=KA,KB)
 406       WRITE(IOUT,602)
 407     2 N=0
 408       DO 3  I=1,NR
 409       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 410       N=N+1
 411       IF(N.LT.10) GO TO 3
 412       WRITE(IOUT,602)
 413       N=0
 414     3 CONTINUE
 415     4 IF (KB.EQ.NC) RETURN
 416       KA=KC+1
 417       KC=KC+6
 418       GO TO 1
 419   600 FORMAT (// 9H ROOT NO.,I4,9I11)
 420   601 FORMAT (/5X,10(1PE11.4))
 421   602 FORMAT (2H  )
 422   603 FORMAT (I5,10F11.5)
 423   604 FORMAT (1H1)
 424       end subroutine EIGOUT
 425 !-----------------------------------------------------------------------------
 426       subroutine MATOUT(NC,NR,LM2,LM3,A)
 427
 428 !      implicit real*8 (a-h,o-z)
 429 !      include 'DIMENSIONS'
 430 !      include 'COMMON.IOUNITS'
 431       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 432       real(kind=8) :: A(LM2,LM3)
 433       KA=1
 434       KC=6
 435     1 KB=MIN0(KC,NC)
 436       WRITE(IOUT,600) (I,I=KA,KB)
 437       WRITE(IOUT,602)
 438     2 N=0
 439       DO 3  I=1,NR
 440       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 441       N=N+1
 442       IF(N.LT.10) GO TO 3
 443       WRITE(IOUT,602)
 444       N=0
 445     3 CONTINUE
 446     4 IF (KB.EQ.NC) RETURN
 447       KA=KC+1
 448       KC=KC+6
 449       GO TO 1
 450   600 FORMAT (//5x,9I11)
 451   602 FORMAT (2H  )
 452   603 FORMAT (I5,10F11.3)
 453   604 FORMAT (1H1)
 454       end subroutine MATOUT
 455 !-----------------------------------------------------------------------------
 456       subroutine MATOUT1(NC,NR,LM2,LM3,A)
 457
 458 !      implicit real*8 (a-h,o-z)
 459 !      include 'DIMENSIONS'
 460 !      include 'COMMON.IOUNITS'
 461       integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N
 462       real(kind=8) :: A(LM2,LM3)
 463       KA=1
 464       KC=21
 465     1 KB=MIN0(KC,NC)
 466       WRITE(IOUT,600) (I,I=KA,KB)
 467       WRITE(IOUT,602)
 468     2 N=0
 469       DO 3  I=1,NR
 470       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 471       N=N+1
 472       IF(N.LT.3) GO TO 3
 473       WRITE(IOUT,602)
 474       N=0
 475     3 CONTINUE
 476     4 IF (KB.EQ.NC) RETURN
 477       KA=KC+1
 478       KC=KC+21
 479       GO TO 1
 480   600 FORMAT (//5x,7(3I5,2x))
 481   602 FORMAT (2H  )
 482   603 FORMAT (I5,7(3F5.1,2x))
 483   604 FORMAT (1H1)
 484       end subroutine MATOUT1
 485 !-----------------------------------------------------------------------------
 486       subroutine MATOUT2(NC,NR,LM2,LM3,A)
 487
 488 !      implicit real*8 (a-h,o-z)
 489 !      include 'DIMENSIONS'
 490 !      include 'COMMON.IOUNITS'
 491       integer :: I,J,KA,KC,KB,N
 492       integer :: LM2,LM3,NC,NR
 493       real(kind=8) :: A(LM2,LM3)
 494       KA=1
 495       KC=12
 496     1 KB=MIN0(KC,NC)
 497       WRITE(IOUT,600) (I,I=KA,KB)
 498       WRITE(IOUT,602)
 499     2 N=0
 500       DO 3  I=1,NR
 501       WRITE(IOUT,603) I,(A(I,J),J=KA,KB)
 502       N=N+1
 503       IF(N.LT.3) GO TO 3
 504       WRITE(IOUT,602)
 505       N=0
 506     3 CONTINUE
 507     4 IF (KB.EQ.NC) RETURN
 508       KA=KC+1
 509       KC=KC+12
 510       GO TO 1
 511   600 FORMAT (//5x,4(3I9,2x))
 512   602 FORMAT (2H  )
 513   603 FORMAT (I5,4(3F9.3,2x))
 514   604 FORMAT (1H1)
 515       end subroutine MATOUT2
 516 !-----------------------------------------------------------------------------
 517       subroutine ginv_mult(z,d_a_tmp)
 518
 519       use geometry_data, only: nres
 520       use control_data
 521       use MPI_data
 522 !      implicit real*8 (a-h,o-z)
 523 !      include 'DIMENSIONS'
 524 #ifdef MPI
 525       include 'mpif.h'
 526       integer :: ierr,ierror
 527 #endif
 528 !      include 'COMMON.SETUP'
 529 !      include 'COMMON.TIME1'
 530 !      include 'COMMON.MD'
 531       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
 532       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
 533       real(kind=8) :: time00,time01
 534       integer :: i,j,k,ind
 535 #ifdef MPI
 536       if (nfgtasks.gt.1) then
 537         if (fg_rank.eq.0) then
 538 ! The matching BROADCAST for fg processors is called in ERGASTULUM
 539           time00=MPI_Wtime()
 540           call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR)
 541           time_Bcast=time_Bcast+MPI_Wtime()-time00
 542 !          print *,"Processor",myrank," BROADCAST iorder in GINV_MULT"
 543         endif
 544 !        write (2,*) "time00",time00
 545 !        write (2,*) "Before Scatterv"
 546 !        call flush(2)
 547 !        write (2,*) "Whole z (for FG master)"
 548 !        do i=1,dimen
 549 !          write (2,*) i,z(i)
 550 !        enddo
 551 !        call MPI_Barrier(FG_COMM,IERROR)
 552         time00=MPI_Wtime()
 553 !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult"
 554         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
 555           MPI_DOUBLE_PRECISION,&
 556           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 557 !        write (2,*) "My chunk of z"
 558 !        do i=1,3*my_ng_count
 559 !          write (2,*) i,z(i)
 560 !        enddo
 561 !        write (2,*) "After SCATTERV"
 562 !        call flush(2)
 563 !        write (2,*) "MPI_Wtime",MPI_Wtime()
 564         time_scatter=time_scatter+MPI_Wtime()-time00
 565 #ifdef TIMING
 566         time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00
 567 #endif
 568 !        write (2,*) "time_scatter",time_scatter
 569 !        write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count",
 570 !     &    my_ng_count
 571 !        call flush(2)
 572         time01=MPI_Wtime()
 573         do k=0,2
 574           do i=1,dimen
 575             ind=(i-1)*3+k+1
 576             temp(ind)=0.0d0
 577             do j=1,my_ng_count
 578 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,
 579 !     &         Ginv(i,j),z((j-1)*3+k+1),
 580 !     &          Ginv(i,j)*z((j-1)*3+k+1)
 581 !              temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1)
 582               temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1)
 583             enddo
 584           enddo
 585         enddo
 586         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 587 !        write (2,*) "Before REDUCE"
 588 !        call flush(2)
 589 !        write (2,*) "z before reduce"
 590 !        do i=1,dimen
 591 !          write (2,*) i,temp(i)
 592 !        enddo
 593         time00=MPI_Wtime()
 594         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
 595             MPI_SUM,king,FG_COMM,IERR)
 596         time_reduce=time_reduce+MPI_Wtime()-time00
 597 !        write (2,*) "After REDUCE"
 598 !        call flush(2)
 599       else
 600 #endif
 601 #ifdef TIMING
 602         time01=MPI_Wtime()
 603 #endif
 604         do k=0,2
 605           do i=1,dimen
 606             ind=(i-1)*3+k+1
 607             d_a_tmp(ind)=0.0d0
 608             do j=1,dimen
 609 !              write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
 610 !              call flush(2)
 611 !     &         Ginv(i,j),z((j-1)*3+k+1),
 612 !     &          Ginv(i,j)*z((j-1)*3+k+1)
 613               d_a_tmp(ind)=d_a_tmp(ind) &
 614                                +Ginv(j,i)*z((j-1)*3+k+1)
 615 !              d_a_tmp(ind)=d_a_tmp(ind)
 616 !     &                         +Ginv(i,j)*z((j-1)*3+k+1)
 617             enddo
 618           enddo
 619         enddo
 620 #ifdef TIMING
 621         time_ginvmult=time_ginvmult+MPI_Wtime()-time01
 622 #endif
 623 #ifdef MPI
 624       endif
 625 #endif
 626       return
 627       end subroutine ginv_mult
 628 !-----------------------------------------------------------------------------
 629 #ifdef GINV_MULT
 630       subroutine ginv_mult_test(z,d_a_tmp)
 631
 632 !      include 'DIMENSIONS'
 633 !el      integer :: dimen
 634 !      include 'COMMON.MD'
 635       real(kind=8),dimension(dimen) :: z,d_a_tmp
 636       real(kind=8),dimension(dimen/3) :: ztmp,dtmp
 637       integer :: i,j,k,ind
 638 !      do i=1,dimen
 639 !        d_a_tmp(i)=0.0d0
 640 !        do j=1,dimen
 641 !          d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j)
 642 !        enddo
 643 !      enddo
 644 !
 645 !      return
 646
 647 !ibm* unroll(3)
 648       do k=0,2
 649        do j=1,dimen/3
 650         ztmp(j)=z((j-1)*3+k+1)
 651        enddo
 652
 653        call alignx(16,ztmp(1))
 654        call alignx(16,dtmp(1))
 655        call alignx(16,Ginv(1,1))
 656
 657        do i=1,dimen/3
 658         dtmp(i)=0.0d0
 659         do j=1,dimen/3
 660            dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j)
 661         enddo
 662        enddo
 663        do i=1,dimen/3
 664         ind=(i-1)*3+k+1
 665         d_a_tmp(ind)=dtmp(i)
 666        enddo
 667       enddo
 668       return
 669       end subroutine ginv_mult_test
 670 #endif
 671 !-----------------------------------------------------------------------------
 672       subroutine fricmat_mult(z,d_a_tmp)
 673
 674       use geometry_data, only: nres
 675       use control_data
 676       use MPI_data
 677 !      include 'DIMENSIONS'
 678 #ifdef MPI
 679       include 'mpif.h'
 680       integer :: IERROR,ierr
 681 #endif
 682 !      include 'COMMON.MD'
 683 !      include 'COMMON.IOUNITS'
 684 !      include 'COMMON.SETUP'
 685 !      include 'COMMON.TIME1'
 686 !#ifndef LANG0
 687 !      include 'COMMON.LANGEVIN'
 688 !#else
 689 !      include 'COMMON.LANGEVIN.lang0'
 690 !#endif
 691       real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp
 692       real(kind=8),dimension(6*nres) :: temp    !(maxres6) maxres6=6*maxres
 693       real(kind=8) :: time00,time01
 694       integer :: i,j,k,ind,nres2
 695       nres2=2*nres
 696 !el      if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
 697
 698 #ifdef MPI
 699       if (nfgtasks.gt.1) then
 700         if (fg_rank.eq.0) then
 701 ! The matching BROADCAST for fg processors is called in ERGASTULUM
 702           time00=MPI_Wtime()
 703           call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR)
 704           time_Bcast=time_Bcast+MPI_Wtime()-time00
 705 !          print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT"
 706         endif
 707 !        call MPI_Barrier(FG_COMM,IERROR)
 708         time00=MPI_Wtime()
 709         call MPI_Scatterv(z,ng_counts(0),ng_start(0),&
 710           MPI_DOUBLE_PRECISION,&
 711           z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
 712 !        write (2,*) "My chunk of z"
 713 !        do i=1,3*my_ng_count
 714 !          write (2,*) i,z(i)
 715 !        enddo
 716         time_scatter=time_scatter+MPI_Wtime()-time00
 717 #ifdef TIMING
 718         time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00
 719 #endif
 720         time01=MPI_Wtime()
 721         do k=0,2
 722           do i=1,dimen
 723             ind=(i-1)*3+k+1
 724             temp(ind)=0.0d0
 725             do j=1,my_ng_count
 726               temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1)
 727             enddo
 728           enddo
 729         enddo
 730         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 731 !        write (2,*) "Before REDUCE"
 732 !        write (2,*) "d_a_tmp before reduce"
 733 !        do i=1,dimen3
 734 !          write (2,*) i,temp(i)
 735 !        enddo
 736 !        call flush(2)
 737         time00=MPI_Wtime()
 738         call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,&
 739             MPI_SUM,king,FG_COMM,IERR)
 740         time_reduce=time_reduce+MPI_Wtime()-time00
 741 !        write (2,*) "After REDUCE"
 742 !        call flush(2)
 743       else
 744 #endif
 745 #ifdef TIMING
 746         time01=MPI_Wtime()
 747 #endif
 748         do k=0,2
 749          do i=1,dimen
 750           ind=(i-1)*3+k+1
 751           d_a_tmp(ind)=0.0d0
 752           do j=1,dimen
 753              d_a_tmp(ind)=d_a_tmp(ind) &
 754                                  -fricmat(j,i)*z((j-1)*3+k+1)
 755           enddo
 756          enddo
 757         enddo
 758 #ifdef TIMING
 759         time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
 760 #endif
 761 #ifdef MPI
 762       endif
 763 #endif
 764 !      write (iout,*) "Vector d_a"
 765 !      do i=1,dimen3
 766 !        write (2,*) i,d_a_tmp(i)
 767 !      enddo
 768       return
 769       end subroutine fricmat_mult
 770 !-----------------------------------------------------------------------------
 771 !-----------------------------------------------------------------------------
 772       end module REMD