+++ /dev/null
- 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
projects / unres4.git / blobdiff