module REMD !----------------------------------------------------------------------------- use io_units use MD_data use REMD_data use muca_md implicit none !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! lagrangian_lesyng.F !----------------------------------------------------------------------------- subroutine lagrangian !------------------------------------------------------------------------- ! This subroutine contains the total lagrangain from which the accelerations ! are obtained. For numerical gradient checking, the derivetive of the ! lagrangian in the velocities and coordinates are calculated seperately !------------------------------------------------------------------------- ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use comm_cipiszcze use energy_data use geometry_data, only: nres use control_data !el, only: mucadyn,lmuca #ifdef MPI include 'mpif.h' real(kind=8) :: time00 #endif ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.MD' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTROL' ! include 'COMMON.MUCA' ! include 'COMMON.TIME1' integer :: i,j,ind,itime real(kind=8) :: zapas(6*nres) !,muca_factor !maxres6=6*maxres logical :: lprn = .false. !el common /cipiszcze/ itime itime = itt_comm #ifdef TIMING time00=MPI_Wtime() #endif do j=1,3 zapas(j)=-gcart(j,0) enddo ind=3 if (lprn) then write (iout,*) "Potential forces backbone" endif do i=nnt,nct-1 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') & i,(-gcart(j,i),j=1,3) do j=1,3 ind=ind+1 zapas(ind)=-gcart(j,i) enddo enddo if (lprn) write (iout,*) "Potential forces sidechain" do i=nnt,nct if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') & i,(-gcart(j,i),j=1,3) do j=1,3 ind=ind+1 zapas(ind)=-gxcart(j,i) enddo endif enddo call ginv_mult(zapas,d_a_work) do j=1,3 d_a(j,0)=d_a_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 d_a(j,i)=d_a_work(ind) enddo enddo do i=nnt,nct if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 ind=ind+1 d_a(j,i+nres)=d_a_work(ind) enddo endif enddo if(lmuca) then imtime=imtime+1 if(mucadyn.gt.0) call muca_update(potE) factor=muca_factor(potE)*t_bath*Rb !d print *,'lmuca ',factor,potE do j=1,3 d_a(j,0)=d_a(j,0)*factor enddo do i=nnt,nct-1 do j=1,3 d_a(j,i)=d_a(j,i)*factor enddo enddo do i=nnt,nct do j=1,3 d_a(j,i+nres)=d_a(j,i+nres)*factor enddo enddo endif if (lprn) then write(iout,*) 'acceleration 3D' write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3) do i=nnt,nct-1 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3) enddo do i=nnt,nct write (iout,'(i3,3f10.5,3x,3f10.5)') & i+nres,(d_a(j,i+nres),j=1,3) enddo endif #ifdef TIMING time_lagrangian=time_lagrangian+MPI_Wtime()-time00 #endif return end subroutine lagrangian !----------------------------------------------------------------------------- subroutine setup_MD_matrices use geometry_data, only: nres,nside use control_data use MPI_data use energy_data use geometry, only:int_bounds use md_calc ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer :: ierror real(kind=8) :: time00 #endif ! include 'COMMON.SETUP' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' logical :: lprn = .false. logical :: osob real(kind=8) :: dtdi real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres !el real(kind=8),dimension(:),allocatable :: Ghalf !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2)) !el real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2) !el real(kind=8),dimension(:,:),allocatable :: Gcopy real(kind=8),dimension(8*6*nres) :: work !(8*maxres6) integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres !el common /przechowalnia/ Gcopy,Ghalf real(kind=8) :: coeff integer :: i,j,ind,ind1,k,ii,jj,m,m1,ii1,iti,nres2,ierr nres2=2*nres if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2) if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2) ! ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv) ! ! Determine the number of degrees of freedom (dimen) and the number of ! sites (dimen1) dimen=(nct-nnt+1)+nside dimen1=(nct-nnt)+(nct-nnt+1) dimen3=dimen*3 #ifdef MPI if (nfgtasks.gt.1) then time00=MPI_Wtime() call MPI_Bcast(5,1,MPI_INTEGER,king,FG_COMM,IERROR) time_Bcast=time_Bcast+MPI_Wtime()-time00 call int_bounds(dimen,igmult_start,igmult_end) igmult_start=igmult_start-1 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,& ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR) my_ng_count=igmult_end-igmult_start call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,& MPI_INTEGER,FG_COMM,IERROR) write (iout,*) 'Processor:',fg_rank,' CG group',kolor,& ' absolute rank',myrank,' igmult_start',igmult_start,& ' igmult_end',igmult_end,' count',my_ng_count write (iout,*) "ng_start",(ng_start(i),i=0,nfgtasks-1) write (iout,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1) call flush(iout) else #endif igmult_start=1 igmult_end=dimen my_ng_count=dimen #ifdef MPI endif #endif ! write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3 ! Zeroing out A and fricmat do i=1,dimen do j=1,dimen A(i,j)=0.0D0 enddo enddo ! Diagonal elements of the dC part of A and the respective friction coefficients ind=1 ind1=0 do i=nnt,nct-1 ind=ind+1 ind1=ind1+1 coeff=0.25d0*IP massvec(ind1)=mp Gmat(ind,ind)=coeff A(ind1,ind)=0.5d0 enddo ! Off-diagonal elements of the dC part of A k=3 do i=1,nct-nnt do j=1,i A(i,j)=1.0d0 enddo enddo ! Diagonal elements of the dX part of A and the respective friction coefficients m=nct-nnt m1=nct-nnt+1 ind=0 ind1=0 msc(ntyp1)=1.0d0 do i=nnt,nct ind=ind+1 ii = ind+m iti=itype(i) massvec(ii)=msc(iabs(iti)) if (iti.ne.10 .and. iti.ne.ntyp1) then ind1=ind1+1 ii1= ind1+m1 A(ii,ii1)=1.0d0 Gmat(ii1,ii1)=ISC(iabs(iti)) endif enddo ! Off-diagonal elements of the dX part of A ind=0 k=nct-nnt do i=nnt,nct iti=itype(i) ind=ind+1 do j=nnt,i ii = ind jj = j-nnt+1 A(k+ii,jj)=1.0d0 enddo enddo if (lprn) then write (iout,*) write (iout,*) "Vector massvec" do i=1,dimen1 write (iout,*) i,massvec(i) enddo write (iout,'(//a)') "A" call matout(dimen,dimen1,nres2,nres2,A) endif ! Calculate the G matrix (store in Gmat) do k=1,dimen do i=1,dimen dtdi=0.0d0 do j=1,dimen1 dtdi=dtdi+A(j,k)*A(j,i)*massvec(j) enddo Gmat(k,i)=Gmat(k,i)+dtdi enddo enddo if (lprn) then write (iout,'(//a)') "Gmat" call matout(dimen,dimen,nres2,nres2,Gmat) endif do i=1,dimen do j=1,dimen Ginv(i,j)=0.0d0 Gcopy(i,j)=Gmat(i,j) enddo Ginv(i,i)=1.0d0 enddo ! Invert the G matrix call MATINVERT(dimen,nres2,Gcopy,Ginv,osob) if (lprn) then write (iout,'(//a)') "Ginv" call matout(dimen,dimen,nres2,nres2,Ginv) endif #ifdef MPI if (nfgtasks.gt.1) then myginv_ng_count=nres2*my_ng_count call MPI_Allgather(nres2*igmult_start,1,MPI_INTEGER,& nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR) call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,& nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR) write (iout,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1) write (iout,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1) call flush(iout) ! call MPI_Scatterv(ginv(1,1),nginv_counts(0), ! & nginv_start(0),MPI_DOUBLE_PRECISION,ginv, ! & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) ! call MPI_Barrier(FG_COMM,IERR) time00=MPI_Wtime() call MPI_Scatterv(ginv(1,1),nginv_counts(0),& nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),& myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) #ifdef TIMING time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00 #endif do i=1,dimen do j=1,2*my_ng_count ginv(j,i)=gcopy(i,j) enddo enddo ! write (iout,*) "Master's chunk of ginv" ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv) endif #endif if (osob) then write (iout,*) "The G matrix is singular." stop endif ! Compute G**(-1/2) and G**(1/2) ind=0 do i=1,dimen do j=1,i ind=ind+1 Ghalf(ind)=Gmat(i,j) enddo enddo call gldiag(nres2,dimen,dimen,Ghalf,work,Geigen,Gvec,& ierr,iwork) if (lprn) then write (iout,'(//a)') & "Eigenvectors and eigenvalues of the G matrix" call eigout(dimen,dimen,nres2,nres2,Gvec,Geigen) endif do i=1,dimen sqreig(i)=dsqrt(Geigen(i)) enddo do i=1,dimen do j=1,dimen Gsqrp(i,j)=0.0d0 Gsqrm(i,j)=0.0d0 Gcopy(i,j)=0.0d0 do k=1,dimen Gsqrp(i,j)=Gsqrp(i,j)+Gvec(i,k)*Gvec(j,k)*sqreig(k) Gsqrm(i,j)=Gsqrm(i,j)+Gvec(i,k)*Gvec(j,k)/sqreig(k) Gcopy(i,j)=Gcopy(i,j)+Gvec(i,k)*Gvec(j,k)*Geigen(k) enddo enddo enddo if (lprn) then write (iout,*) "Comparison of original and restored G" do i=1,dimen do j=1,dimen write (iout,'(2i5,5f10.5)') i,j,Gmat(i,j),Gcopy(i,j),& Gmat(i,j)-Gcopy(i,j),Gsqrp(i,j),Gsqrm(i,j) enddo enddo endif ! deallocate(Gcopy) return end subroutine setup_MD_matrices !----------------------------------------------------------------------------- subroutine EIGOUT(NC,NR,LM2,LM3,A,B) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N real(kind=8) :: A(LM2,LM3),B(LM2) KA=1 KC=6 1 KB=MIN0(KC,NC) WRITE(IOUT,600) (I,I=KA,KB) WRITE(IOUT,601) (B(I),I=KA,KB) WRITE(IOUT,602) 2 N=0 DO 3 I=1,NR WRITE(IOUT,603) I,(A(I,J),J=KA,KB) N=N+1 IF(N.LT.10) GO TO 3 WRITE(IOUT,602) N=0 3 CONTINUE 4 IF (KB.EQ.NC) RETURN KA=KC+1 KC=KC+6 GO TO 1 600 FORMAT (// 9H ROOT NO.,I4,9I11) 601 FORMAT (/5X,10(1PE11.4)) 602 FORMAT (2H ) 603 FORMAT (I5,10F11.5) 604 FORMAT (1H1) end subroutine EIGOUT !----------------------------------------------------------------------------- subroutine MATOUT(NC,NR,LM2,LM3,A) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N real(kind=8) :: A(LM2,LM3) KA=1 KC=6 1 KB=MIN0(KC,NC) WRITE(IOUT,600) (I,I=KA,KB) WRITE(IOUT,602) 2 N=0 DO 3 I=1,NR WRITE(IOUT,603) I,(A(I,J),J=KA,KB) N=N+1 IF(N.LT.10) GO TO 3 WRITE(IOUT,602) N=0 3 CONTINUE 4 IF (KB.EQ.NC) RETURN KA=KC+1 KC=KC+6 GO TO 1 600 FORMAT (//5x,9I11) 602 FORMAT (2H ) 603 FORMAT (I5,10F11.3) 604 FORMAT (1H1) end subroutine MATOUT !----------------------------------------------------------------------------- subroutine MATOUT1(NC,NR,LM2,LM3,A) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: LM2,LM3,NC,NR,KA,KC,KB,I,J,N real(kind=8) :: A(LM2,LM3) KA=1 KC=21 1 KB=MIN0(KC,NC) WRITE(IOUT,600) (I,I=KA,KB) WRITE(IOUT,602) 2 N=0 DO 3 I=1,NR WRITE(IOUT,603) I,(A(I,J),J=KA,KB) N=N+1 IF(N.LT.3) GO TO 3 WRITE(IOUT,602) N=0 3 CONTINUE 4 IF (KB.EQ.NC) RETURN KA=KC+1 KC=KC+21 GO TO 1 600 FORMAT (//5x,7(3I5,2x)) 602 FORMAT (2H ) 603 FORMAT (I5,7(3F5.1,2x)) 604 FORMAT (1H1) end subroutine MATOUT1 !----------------------------------------------------------------------------- subroutine MATOUT2(NC,NR,LM2,LM3,A) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: I,J,KA,KC,KB,N integer :: LM2,LM3,NC,NR real(kind=8) :: A(LM2,LM3) KA=1 KC=12 1 KB=MIN0(KC,NC) WRITE(IOUT,600) (I,I=KA,KB) WRITE(IOUT,602) 2 N=0 DO 3 I=1,NR WRITE(IOUT,603) I,(A(I,J),J=KA,KB) N=N+1 IF(N.LT.3) GO TO 3 WRITE(IOUT,602) N=0 3 CONTINUE 4 IF (KB.EQ.NC) RETURN KA=KC+1 KC=KC+12 GO TO 1 600 FORMAT (//5x,4(3I9,2x)) 602 FORMAT (2H ) 603 FORMAT (I5,4(3F9.3,2x)) 604 FORMAT (1H1) end subroutine MATOUT2 !----------------------------------------------------------------------------- subroutine ginv_mult(z,d_a_tmp) use geometry_data, only: nres use control_data use MPI_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer :: ierr,ierror #endif ! include 'COMMON.SETUP' ! include 'COMMON.TIME1' ! include 'COMMON.MD' real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp real(kind=8),dimension(6*nres) :: temp !(maxres6) maxres6=6*maxres real(kind=8) :: time00,time01 integer :: i,j,k,ind #ifdef MPI if (nfgtasks.gt.1) then if (fg_rank.eq.0) then ! The matching BROADCAST for fg processors is called in ERGASTULUM time00=MPI_Wtime() call MPI_Bcast(4,1,MPI_INTEGER,king,FG_COMM,IERROR) time_Bcast=time_Bcast+MPI_Wtime()-time00 ! print *,"Processor",myrank," BROADCAST iorder in GINV_MULT" endif ! write (2,*) "time00",time00 ! write (2,*) "Before Scatterv" ! call flush(2) ! write (2,*) "Whole z (for FG master)" ! do i=1,dimen ! write (2,*) i,z(i) ! enddo ! call MPI_Barrier(FG_COMM,IERROR) time00=MPI_Wtime() !elwrite(iout,*) "do tej pory jest OK, MPI_Scatterv w ginv_mult" call MPI_Scatterv(z,ng_counts(0),ng_start(0),& MPI_DOUBLE_PRECISION,& z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) ! write (2,*) "My chunk of z" ! do i=1,3*my_ng_count ! write (2,*) i,z(i) ! enddo ! write (2,*) "After SCATTERV" ! call flush(2) ! write (2,*) "MPI_Wtime",MPI_Wtime() time_scatter=time_scatter+MPI_Wtime()-time00 #ifdef TIMING time_scatter_ginvmult=time_scatter_ginvmult+MPI_Wtime()-time00 #endif ! write (2,*) "time_scatter",time_scatter ! write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count", ! & my_ng_count ! call flush(2) time01=MPI_Wtime() do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 temp(ind)=0.0d0 do j=1,my_ng_count ! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1, ! & Ginv(i,j),z((j-1)*3+k+1), ! & Ginv(i,j)*z((j-1)*3+k+1) ! temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1) temp(ind)=temp(ind)+Ginv(j,i)*z1((j-1)*3+k+1) enddo enddo enddo time_ginvmult=time_ginvmult+MPI_Wtime()-time01 ! write (2,*) "Before REDUCE" ! call flush(2) ! write (2,*) "z before reduce" ! do i=1,dimen ! write (2,*) i,temp(i) ! enddo time00=MPI_Wtime() call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,& MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 ! write (2,*) "After REDUCE" ! call flush(2) else #endif #ifdef TIMING time01=MPI_Wtime() #endif do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 d_a_tmp(ind)=0.0d0 do j=1,dimen ! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1 ! call flush(2) ! & Ginv(i,j),z((j-1)*3+k+1), ! & Ginv(i,j)*z((j-1)*3+k+1) d_a_tmp(ind)=d_a_tmp(ind) & +Ginv(j,i)*z((j-1)*3+k+1) ! d_a_tmp(ind)=d_a_tmp(ind) ! & +Ginv(i,j)*z((j-1)*3+k+1) enddo enddo enddo #ifdef TIMING time_ginvmult=time_ginvmult+MPI_Wtime()-time01 #endif #ifdef MPI endif #endif return end subroutine ginv_mult !----------------------------------------------------------------------------- #ifdef GINV_MULT subroutine ginv_mult_test(z,d_a_tmp) ! include 'DIMENSIONS' !el integer :: dimen ! include 'COMMON.MD' real(kind=8),dimension(dimen) :: z,d_a_tmp real(kind=8),dimension(dimen/3) :: ztmp,dtmp integer :: i,j,k,ind ! do i=1,dimen ! d_a_tmp(i)=0.0d0 ! do j=1,dimen ! d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j) ! enddo ! enddo ! ! return !ibm* unroll(3) do k=0,2 do j=1,dimen/3 ztmp(j)=z((j-1)*3+k+1) enddo call alignx(16,ztmp(1)) call alignx(16,dtmp(1)) call alignx(16,Ginv(1,1)) do i=1,dimen/3 dtmp(i)=0.0d0 do j=1,dimen/3 dtmp(i)=dtmp(i)+Ginv(i,j)*ztmp(j) enddo enddo do i=1,dimen/3 ind=(i-1)*3+k+1 d_a_tmp(ind)=dtmp(i) enddo enddo return end subroutine ginv_mult_test #endif !----------------------------------------------------------------------------- subroutine fricmat_mult(z,d_a_tmp) use geometry_data, only: nres use control_data use MPI_data ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer :: IERROR,ierr #endif ! include 'COMMON.MD' ! include 'COMMON.IOUNITS' ! include 'COMMON.SETUP' ! include 'COMMON.TIME1' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif real(kind=8),dimension(dimen3) :: z,z1,d_a_tmp real(kind=8),dimension(6*nres) :: temp !(maxres6) maxres6=6*maxres real(kind=8) :: time00,time01 integer :: i,j,k,ind,nres2 nres2=2*nres !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) #ifdef MPI if (nfgtasks.gt.1) then if (fg_rank.eq.0) then ! The matching BROADCAST for fg processors is called in ERGASTULUM time00=MPI_Wtime() call MPI_Bcast(9,1,MPI_INTEGER,king,FG_COMM,IERROR) time_Bcast=time_Bcast+MPI_Wtime()-time00 ! print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT" endif ! call MPI_Barrier(FG_COMM,IERROR) time00=MPI_Wtime() call MPI_Scatterv(z,ng_counts(0),ng_start(0),& MPI_DOUBLE_PRECISION,& z1,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) ! write (2,*) "My chunk of z" ! do i=1,3*my_ng_count ! write (2,*) i,z(i) ! enddo time_scatter=time_scatter+MPI_Wtime()-time00 #ifdef TIMING time_scatter_fmatmult=time_scatter_fmatmult+MPI_Wtime()-time00 #endif time01=MPI_Wtime() do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 temp(ind)=0.0d0 do j=1,my_ng_count temp(ind)=temp(ind)-fricmat(j,i)*z1((j-1)*3+k+1) enddo enddo enddo time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01 ! write (2,*) "Before REDUCE" ! write (2,*) "d_a_tmp before reduce" ! do i=1,dimen3 ! write (2,*) i,temp(i) ! enddo ! call flush(2) time00=MPI_Wtime() call MPI_Reduce(temp(1),d_a_tmp(1),dimen3,MPI_DOUBLE_PRECISION,& MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 ! write (2,*) "After REDUCE" ! call flush(2) else #endif #ifdef TIMING time01=MPI_Wtime() #endif do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 d_a_tmp(ind)=0.0d0 do j=1,dimen d_a_tmp(ind)=d_a_tmp(ind) & -fricmat(j,i)*z((j-1)*3+k+1) enddo enddo enddo #ifdef TIMING time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01 #endif #ifdef MPI endif #endif ! write (iout,*) "Vector d_a" ! do i=1,dimen3 ! write (2,*) i,d_a_tmp(i) ! enddo return end subroutine fricmat_mult !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module REMD