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 #ifdef FIVEDIAG use energy, only: grad_transform #endif 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,mnum,innt,inct,inct_prev,ichain,n,mark real(kind=8) :: zapas(6*nres) !,muca_factor !maxres6=6*maxres ! the line below might be wrong #ifdef FIVEDIAG real(kind=8) :: rs(2*nres),xsolv(2*nres) !#ifdef CHECK5DSOL real(kind=8) :: rscheck(2*nres),rsold(2*nres) !#endif #endif logical :: lprn = .false. !el common /cipiszcze/ itime itime = itt_comm #ifdef TIMING time00=MPI_Wtime() #endif #ifdef FIVEDIAG call grad_transform d_a=0.0d0 if (lprn) then write (iout,*) "Potential forces backbone" do i=1,nres write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3) enddo write (iout,*) "Potential forces sidechain" do i=nnt,nct ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) & write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3) enddo endif do ichain=1,nchain n=dimen_chain(ichain) innt=iposd_chain(ichain) do j=1,3 ind=1 do i=chain_border(1,ichain),chain_border(2,ichain) mnum=molnum(i) if (itype(i,1).eq.10.or.mnum.ge.3)then rs(ind)=-gcart(j,i)-gxcart(j,i) ind=ind+1 else rs(ind)=-gcart(j,i) rs(ind+1)=-gxcart(j,i) ind=ind+2 end if enddo #ifdef CHECK5DSOL rsold=rs #endif if (lprn) then write(iout,*) "RHS of the 5-diag equations system",& ichain," j",j do i=1,n write(iout,*) i,rs(i) enddo endif call FDISYS (n,DM(innt),DU1(innt),DU2(innt),rs,xsolv) if (lprn) then write (iout,*) "Solution of the 5-diagonal equations system" do i=1,n write (iout,'(i5,f10.5)') i,xsolv(i) enddo endif #ifdef CHECK5DSOL ! Check the solution call fivediagmult(n,DMorig(innt),DU1orig(innt),DU2orig(innt),& xsolv,rscheck) do i=1,n write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),& "ratio",rscheck(i)/rsold(i) enddo ! end check #endif #undef CHECK5DSOL ind=1 do i=chain_border(1,ichain),chain_border(2,ichain) mnum=molnum(i) if (itype(i,1).eq.10 .or.mnum.ge.3) then d_a(j,i)=xsolv(ind) ind=ind+1 else d_a(j,i)=xsolv(ind) d_a(j,i+nres)=xsolv(ind+1) ind=ind+2 end if enddo enddo ! j enddo ! ichain if (lprn) then write (iout,*) "Acceleration in CA and SC oordinates" do i=1,nres write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3) enddo do i=nnt,nct write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3) enddo endif !C Conevert d_a to virtual-bon-vector basis #define WLOS #ifdef WLOS !c write (iout,*) "WLOS" if (nnt.eq.1) then d_a(:,0)=d_a(:,1) endif do i=1,nres mnum=molnum(i) if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum) .or.mnum.ge.3) then do j=1,3 d_a(j,i)=d_a(j,i+1)-d_a(j,i) enddo else do j=1,3 d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i) d_a(j,i)=d_a(j,i+1)-d_a(j,i) enddo end if enddo d_a(:,nres)=0.0d0 d_a(:,nct)=0.0d0 d_a(:,2*nres)=0.0d0 !c d_a(:,0)=d_a(:,1) !c d_a(:,1)=0.0d0 !c write (iout,*) "Shifting accelerations" if (nnt.gt.1) then d_a(:,0)=d_a(:,1) d_a(:,1)=0.0d0 endif #define CHUJ #ifdef CHUJ do ichain=2,nchain !TEST 27.06.2023 godz 16.00 if (molnum(chain_border1(1,ichain)+1).eq.5) cycle !c write (iout,*) "ichain",chain_border1(1,ichain)-1, !c & chain_border1(1,ichain) d_a(:,chain_border1(1,ichain)-1)=d_a(:,chain_border1(1,ichain)) d_a(:,chain_border1(1,ichain))=0.0d0 enddo !c write (iout,*) "Adding accelerations" do ichain=2,nchain if (molnum(chain_border1(1,ichain)+1).eq.5) cycle !c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1, !c & chain_border(2,ichain-1) d_a(:,chain_border1(1,ichain)-1)=& d_a(:,chain_border1(1,ichain)-1)+d_a(:,chain_border(2,ichain-1)) d_a(:,chain_border(2,ichain-1))=0.0d0 enddo #endif #else inct_prev=0 do j=1,3 aaux(j)=0.0d0 enddo do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) do j=1,3 d_a(j,inct_prev)=d_a(j,innt)-aaux(j) enddo inct_prev=inct+1 do i=innt,inct mnum=molnum(i) if ((itype(i).ne.10).and.(mnum.le.3)) then do j=1,3 d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i) enddo endif enddo do j=1,3 aaux(j)=d_a(j,inct) enddo do i=innt,inct do j=1,3 d_a(j,i)=d_a(j,i+1)-d_a(j,i) enddo enddo enddo #endif if (lprn) then write(iout,*) 'acceleration 3D FIVEDIAG in dC and dX' do i=0,nres 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,(d_a(j,i+nres),j=1,3) enddo endif #else ! Old procedure 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 mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.lt.4) then if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') & i,(-gxcart(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 mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.lt.4) then do j=1,3 ind=ind+1 d_a(j,i+nres)=d_a_work(ind) enddo endif enddo #endif 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 ! real(kind=8) ,allocatable, dimension(:) :: DDM,DDU1,DDU2 #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 #ifndef FIVEDIAG real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult #endif real(kind=8) :: dtdi real(kind=8),dimension(:),allocatable :: massvec,sqreig !(maxres2) maxres2=2*maxres real(kind=8) :: relfeh,eps1,eps2 !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(:),allocatable :: work !(8*maxres6) integer,dimension(:),allocatable :: iwork !(maxres6) maxres6=6*maxres ! common /jakistam/ iwork,work,massvec,sqreig !el common /przechowalnia/ Gcopy,Ghalf real(kind=8) :: coeff,mscab integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark real(kind=8) :: ip4 integer :: iz,mnum,ichain,n,dimenp,innt,inct if(.not.allocated(massvec)) then allocate(massvec(2*nres),sqreig(2*nres)) allocate(work(8*6*nres)) allocate(iwork(6*nres)) endif print *,"just entered" relfeh=1.0d-14 nres2=2*nres print *,"FIVEDIAG" write (iout,*) "before FIVEDIAG" #ifdef FIVEDIAG if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2) #endif #ifndef FIVEDIAG if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2) write (iout,*) "ALLOCATE" print *,"ALLOCATE" 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) if(.not.allocated(Bmat)) allocate(Bmat(nres2,nres2)) if(.not.allocated(matmult)) allocate(matmult(nres2,nres2)) #endif ! ! 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 ! write (iout,*) "nnt",nnt," nct",nct," nside",nside #ifdef FIVEDIAG #ifdef CRYST_BOND ip4=ip/4 DM(1)=mp/4+ip4 if (iabs(itype(nnt)).eq.10) then DM(1)=DM(1)+msc(10) ind=2 nind=1 else DM(1)=DM(1)+isc(iabs(itype(nnt))) DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt))) ind=3 nind=2 endif do i=nnt+1,nct-1 mnum=molnum(i) ! if (iabs(itype(i,1)).eq.ntyp1) cycle DM(ind)=2*ip4+mp/2 if (iabs(itype(i,1)).eq.10 .or. & iabs(itype(i,mnum)).eq.ntyp1_molec(mnum) .or. mnum.ge.4) then if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10) ind=ind+1 else DM(ind)=DM(ind)+isc(iabs(itype(i,1))) DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1))) ind=ind+2 endif enddo if (nct.gt.nnt) then DM(ind)=ip4+mp/4 if (iabs(itype(nct)).eq.10) then DM(ind)=DM(ind)+msc(10) nind=ind else DM(ind)=DM(ind)+isc(iabs(itype(nct))) DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct))) nind=ind+1 endif endif ind=1 do i=nnt,nct mnum=molnum(i) if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 & .and.mnum.lt.4) then DU1(ind)=-isc(iabs(itype(i,1))) DU1(ind+1)=0.0d0 ind=ind+2 else DU1(ind)=mp/4-ip4 ind=ind+1 endif enddo ind=1 do i=nnt,nct-1 mnum=molnum(i) ! if (iabs(itype(i,1)).eq.ntyp1) cycle write (iout,*) "i",i," itype",itype(i,1),ntyp1 if (iabs(itype(i,1)).ne.10 .and. & iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.lt.4) then DU2(ind)=mp/4-ip4 DU2(ind+1)=0.0d0 ind=ind+2 else DU2(ind)=0.0d0 DU2(ind+1)=0.0d0 ind=ind+1 endif enddo #else dimen=0 dimen1=0 dimenp=0 do ichain=1,nchain dimen=dimen+chain_length(ichain) dimen1=dimen1+2*chain_length(ichain)-1 dimenp=dimenp+chain_length(ichain)-1 enddo write (iout,*) "Number of Calphas",dimen write (iout,*) "Number of sidechains",nside write (iout,*) "Number of peptide groups",dimenp dimen=dimen+nside ! number of centers dimen3=3*dimen ! degrees of freedom write (iout,*) "Number of centers",dimen write (iout,*) "Degrees of freedom:",dimen3 ! ip4=ip/4 ind=1 do ichain=1,nchain iposd_chain(ichain)=ind innt=chain_border(1,ichain) mnum=molnum(innt) inct=chain_border(2,ichain) if (mnum.eq.5) mp(mnum)=0.0 if (mnum.eq.5) ip(mnum)=0.0 ! if (mnum.eq.5) mp(mnum)=msc(itype(innt,mnum),mnum) DM(ind)=mp(mnum)/4+ip(mnum)/4 if (iabs(itype(innt,1)).eq.10.or.molnum(innt).gt.2) then DM(ind)=DM(ind)+msc(itype(innt,mnum),mnum) ind=ind+1 nind=1 else DM(ind)=DM(ind)+isc(iabs(itype(innt,mnum)),mnum) DM(ind+1)=msc(iabs(itype(innt,mnum)),mnum)+isc(iabs(itype(innt,mnum)),mnum) ind=ind+2 nind=2 endif write (iout,*) "ind",ind," nind",nind do i=innt+1,inct-1 mnum=molnum(i) ! if (iabs(itype(i)).eq.ntyp1) cycle ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) ! if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum) if (mnum.eq.5) mp(mnum)=0.0 if (mnum.eq.5) ip(mnum)=0.0 ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) DM(ind)=2*ip(mnum)/4+mp(mnum)/2 if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2) then ! if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2)& ! DM(ind)=DM(ind)+msc(itype(i,molnum(i)),mnum) DM(ind)=DM(ind)+msc(itype(i,mnum),mnum) ind=ind+1 nind=nind+1 else DM(ind)=DM(ind)+isc(iabs(itype(i,mnum)),mnum) DM(ind+1)=msc(iabs(itype(i,mnum)),mnum)+isc(iabs(itype(i,mnum)),mnum) ind=ind+2 nind=nind+2 endif write (iout,*) "i",i," ind",ind," nind",nind enddo if (inct.gt.innt) then ! DM(ind)=ip4+mp(molnum(inct))/4 mnum=molnum(inct) if (mnum.eq.5) mp(mnum)=0.0 ! if (mnum.eq.5) mp(mnum)=msc(itype(inct,molnum(inct)),molnum(inct)) DM(ind)=mp(mnum)/4+ip(mnum)/4 if (iabs(itype(inct,mnum)).eq.10.or.molnum(inct).gt.2) then DM(ind)=DM(ind)+msc(itype(inct,molnum(inct)),molnum(inct)) ind=ind+1 nind=nind+1 else mnum=molnum(inct) DM(ind)=DM(ind)+isc(iabs(itype(inct,mnum)),mnum) DM(ind+1)=msc(iabs(itype(inct,mnum)),mnum)+isc(iabs(itype(inct,mnum)),mnum) ind=ind+2 nind=nind+2 endif endif write (iout,*) "ind",ind," nind",nind dimen_chain(ichain)=nind enddo do ichain=1,nchain ind=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=innt,inct mnum=molnum(i) if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then DU1(ind)=-isc(iabs(itype(i,mnum)),mnum) DU1(ind+1)=0.0d0 ind=ind+2 else ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) ! if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum) if (mnum.eq.5) mp(mnum)=0.0 DU1(ind)=mp(mnum)/4-ip(mnum)/4 ind=ind+1 endif enddo enddo do ichain=1,nchain ind=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=innt,inct-1 mnum=molnum(i) ! if (iabs(itype(i)).eq.ntyp1) cycle !c write (iout,*) "i",i," itype",itype(i),ntyp1 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) if (mnum.eq.5) mp(mnum)=0.0 DU2(ind)=mp(mnum)/4-ip(mnum)/4 DU2(ind+1)=0.0d0 ind=ind+2 else DU2(ind)=0.0d0 DU2(ind+1)=0.0d0 ind=ind+1 endif enddo enddo DMorig=DM DU1orig=DU1 DU2orig=DU2 gmatout=.true. if (gmatout) then write (iout,*)"The upper part of the five-diagonal inertia matrix" endif do ichain=1,nchain if (gmatout) write (iout,'(a,i5)') 'Chain',ichain n=dimen_chain(ichain) innt=iposd_chain(ichain) inct=iposd_chain(ichain)+dimen_chain(ichain)-1 if (gmatout) then do i=innt,inct if (i.lt.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i) else if (i.eq.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i) else write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i) endif enddo endif call FDISYP (n, DM(innt:inct), DU1(innt:inct-1),& DU2(innt:inct-1), MARK) if (mark.eq.-1) then write(iout,*)& "ERROR: the inertia matrix is not positive definite for chain",& ichain #ifdef MPI call MPI_Finalize(ierr) #endif stop else if (mark.eq.0) then write (iout,*)& "ERROR: the inertia matrix is singular for chain",ichain #ifdef MPI call MPI_Finalize(ierr) #endif else if (mark.eq.1) then if (gmatout) then write (iout,*) "The transformed five-diagonal inertia matrix" write (iout,'(a,i5)') 'Chain',ichain do i=innt,inct if (i.lt.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i),DU2(i) else if (i.eq.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i),DU1(i) else write (iout,'(2i3,3f10.5)') i,i-innt+1,DM(i) endif enddo endif endif enddo ! Diagonalization of the pentadiagonal matrix #ifdef TIMING time00=MPI_Wtime() #endif #endif !CRYSTBOND #else dimen=(nct-nnt+1)+nside dimen1=(nct-nnt)+(nct-nnt+1) dimen3=dimen*3 write (iout,*) "nnt",nnt," nct",nct," nside",nside ! Old Gmatrix #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 ! print *,"TUTUTUT?!",nnt,nct-1 do i=nnt,nct-1 mnum=molnum(i) ind=ind+1 ind1=ind1+1 coeff=0.25d0*IP(mnum) print *,"TU",i,itype(i,mnum),mnum if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) print *,"TU2",i,coeff,mp(mnum) massvec(ind1)=mp(mnum) Gmat(ind,ind)=coeff ! print *,"i",mp(mnum) 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 do i=1,2 msc(ntyp1_molec(i),i)=1.0d0 enddo msc(ntyp1_molec(4),4)=1.0d0 ! print *,"TU3",ntyp1_molec(4)-1 msc(ntyp1_molec(4)-1,4)=1.0d0 do i=nnt,nct ind=ind+1 ii = ind+m mnum=molnum(i) iti=itype(i,mnum) ! if (mnum.eq.5) then ! mscab=0.0 ! else mscab=msc(iabs(iti),mnum) ! endif massvec(ii)=mscab if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.lt.4) then ind1=ind1+1 ii1= ind1+m1 A(ii,ii1)=1.0d0 Gmat(ii1,ii1)=ISC(iabs(iti),mnum) endif enddo ! Off-diagonal elements of the dX part of A ind=0 k=nct-nnt do i=nnt,nct mnum=molnum(i) iti=itype(i,mnum) 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 Bmat=0.0d0 Bmat(1,1)=1.0d0 j=2 do i=nnt,nct mnum=molnum(i) if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))& .or.(mnum.ge.4)) then Bmat(i-nnt+2,j-1)=-1 Bmat(i-nnt+2,j)=1 j=j+1 else if (i.lt.nct) then Bmat(i-nnt+2,j-1)=-1 Bmat(i-nnt+2,j+1)=1 Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1 Bmat(i-nnt+1+(nct-nnt)+1,j)=1 j=j+2 else Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1 Bmat(i-nnt+1+(nct-nnt)+1,j)=1 endif endif enddo write (iout,*) "j",j," dimen",dimen if (lprn) then write (iout,'(//a)') "Bmat" call matout(dimen,dimen,nres2,nres2,Bmat) endif Gcopy=0.0d0 write(iout,*) "before Gcopy",dimen,nres*2 #ifdef TESTFIVEDIAG do i=1,dimen do j=1,dimen ! write (iout,*) "myij",i,j do k=1,dimen do l=1,dimen ! write(iout,*) "i,j,k,l",i,j,k,l Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j) enddo enddo enddo enddo #endif if (lprn) then write (iout,'(//a)') "Gmat-transformed" call matout(dimen,dimen,nres2,nres2,Gcopy) 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." write (iout,'(//a)') "Gmat-transformed" call matout(dimen,dimen,nres2,nres2,Gcopy) 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 if (allocated(Gcopy)) deallocate(Gcopy) #endif !write(iout,*) "end setup_MD_matr" 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 #ifdef FIVEDIAG subroutine fivediagmult(n,DM,DU1,DU2,x,y) integer n double precision DM(n),DU1(n),DU2(n),x(n),y(n) integer i y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3) y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4) do i=3,n-2 y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i) & +DU1(i)*x(i+1)+DU2(i)*x(i+2) enddo y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1) & +DU1(n-1)*x(n) y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n) return end subroutine !c--------------------------------------------------------------------------- subroutine fivediaginv_mult(ndim,forces,d_a_vec) use energy_data, only:nchain,chain_border,nct,nnt,molnum,& chain_border1,itype use geometry_data, only: nside integer ndim real(kind=8),dimension(:),allocatable ::forces,d_a_vec real(kind=8),dimension(:),allocatable :: xsolv,rs real(kind=8),dimension(:,:),allocatable :: accel integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev,mnum accel=0.0d0 !#define DEBUG if (.not.allocated(forces)) allocate(forces(6*nres)) if (.not.allocated(d_a_vec)) allocate(d_a_vec(6*nres)) if (.not.allocated(xsolv)) allocate(xsolv(3*ndim)) if (.not.allocated(rs)) allocate(rs(3*ndim)) if (.not.allocated(accel)) allocate(accel(3,0:2*nres)) #ifdef DEBUG do i=1,6*nres write(iout,*) "received forces",i,forces(i) enddo #endif do j=1,3 !Compute accelerations in Calpha and SC do ichain=1,nchain n=dimen_chain(ichain) iposc=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=iposc,iposc+n-1 ! write(iout,*) "index",3*(i-1)+j,forces(3*(i-1)+j) rs(i-iposc+1)=forces(3*(i-1)+j) enddo #ifdef DEBUG write (iout,*) "j",j," chain",ichain,"n",n write (iout,*) "innt",innt,inct,iposc write (iout,*) "rs" do i=1,n write (iout,'(i5,f10.5)') i,rs(i) enddo #endif call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv) #ifdef DEBUG write (iout,*) "xsolv" write (iout,'(f10.5)') (xsolv(i),i=1,n) #endif ind=1 do i=innt,inct mnum=molnum(i) if (itype(i,1).eq.10.or.mnum.gt.2)then accel(j,i)=xsolv(ind) ind=ind+1 else accel(j,i)=xsolv(ind) accel(j,i+nres)=xsolv(ind+1) ind=ind+2 end if enddo enddo enddo !C Convert d_a to virtual-bon-vector basis #ifdef DEBUG write (iout,*) "accel in CA-SC basis" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),& (accel(j,i+nres),j=1,3) enddo write (iout,*) "nnt",nnt #endif if (nnt.eq.1) then accel(:,0)=accel(:,1) endif do i=1,nres mnum=molnum(i) if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& .or.mnum.ge.3) then do j=1,3 accel(j,i)=accel(j,i+1)-accel(j,i) enddo else do j=1,3 accel(j,i+nres)=accel(j,i+nres)-accel(j,i) accel(j,i)=accel(j,i+1)-accel(j,i) enddo end if enddo accel(:,nres)=0.0d0 accel(:,nct)=0.0d0 accel(:,2*nres)=0.0d0 if (nnt.gt.1) then accel(:,0)=accel(:,1) accel(:,1)=0.0d0 endif do ichain=2,nchain if (molnum(chain_border1(1,ichain)+1).eq.5) cycle accel(:,chain_border1(1,ichain)-1)= & accel(:,chain_border1(1,ichain)) accel(:,chain_border1(1,ichain))=0.0d0 enddo do ichain=2,nchain if (molnum(chain_border1(1,ichain)+1).eq.5) cycle accel(:,chain_border1(1,ichain)-1)= & accel(:,chain_border1(1,ichain)-1) & +accel(:,chain_border(2,ichain-1)) accel(:,chain_border(2,ichain-1))=0.0d0 enddo #ifdef DEBUG write (iout,*) "accel in fivediaginv_mult: 1" do i=0,2*nres write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3) enddo #endif do j=1,3 d_a_vec(j)=accel(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_a_vec(ind+j)=accel(j,i) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.le.2) then do j=1,3 d_a_vec(ind+j)=accel(j,i+nres) enddo ind=ind+3 endif enddo #ifdef DEBUG write (iout,*) "d_a_vec" write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside)) #endif return end subroutine #undef DEBUG #else !----------------------------------------------------------------------------- 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 !----------------------------------------------------------------------------- #endif !----------------------------------------------------------------------------- end module REMD