subroutine lagrangian c------------------------------------------------------------------------- c This subroutine contains the total lagrangain from which the accelerations c are obtained. For numerical gradient checking, the derivetive of the c lagrangian in the velocities and coordinates are calculated seperately c------------------------------------------------------------------------- implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer time00 #endif include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif #ifdef LANG0 #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #else include 'COMMON.LANGEVIN' #endif include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.MUCA' include 'COMMON.TIME1' integer i,j,ind double precision zapas(MAXRES6),muca_factor logical lprn /.false./ integer itime common /cipiszcze/ itime #ifdef FIVEDIAG double precision rs(maxres2_chain),xsolv(maxres2_chain),ip4 double precision aaux(3) integer nind,innt,inct,inct_prev,ichain,n,mark #ifdef CHECK5DSOL double precision rscheck(maxres2_chain),rsold(maxres2_chain) #endif #endif #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) if (itype(i).eq.10)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, chain",ichain," j",j do i=1,n write(iout,'(i5,f10.5)') 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) if (itype(i).eq.10) 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 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) 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 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 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 if (itype(i).ne.10) 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 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 cd 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 #endif #ifdef TIMING time_lagrangian=time_lagrangian+MPI_Wtime()-time00 #endif return end c------------------------------------------------------------------ subroutine setup_MD_matrices implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierror,ierr double precision time00 #endif include 'COMMON.CONTROL' include 'COMMON.SETUP' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #endif include 'COMMON.IOUNITS' include 'COMMON.TIME1' integer i,j,k,m,m1,ind,ind1,ii,iti,ii1,jj double precision coeff logical lprn /.false./ logical osob double precision dtdi,massvec(maxres2) #ifdef FIVEDIAG integer ichain,innt,inct,nind,mark,n double precision ip4 #else double precision Gcopy(maxres2,maxres2),Ghalf(mmaxres2), & sqreig(maxres2) double precision work(8*maxres6) integer iwork(maxres6) common /przechowalnia/ Gcopy,Ghalf #endif c c Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the c inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv) c c Determine the number of degrees of freedom (dimen) and the number of c sites (dimen1) #ifdef FIVEDIAG dimen=0 dimen1=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) inct=chain_border(2,ichain) DM(ind)=mp/4+ip4 if (iabs(itype(innt)).eq.10) then DM(ind)=DM(ind)+msc(10) ind=ind+1 nind=1 else DM(ind)=DM(ind)+isc(iabs(itype(innt))) DM(ind+1)=msc(iabs(itype(innt)))+isc(iabs(itype(innt))) ind=ind+2 nind=2 endif write (iout,*) "ind",ind," nind",nind do i=innt+1,inct-1 ! if (iabs(itype(i)).eq.ntyp1) cycle DM(ind)=2*ip4+mp/2 if (iabs(itype(i)).eq.10) then if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10) ind=ind+1 nind=nind+1 else DM(ind)=DM(ind)+isc(iabs(itype(i))) DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i))) 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/4 if (iabs(itype(inct)).eq.10) then DM(ind)=DM(ind)+msc(10) ind=ind+1 nind=nind+1 else DM(ind)=DM(ind)+isc(iabs(itype(inct))) DM(ind+1)=msc(iabs(itype(inct)))+isc(iabs(itype(inct))) 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 if (iabs(itype(i)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then DU1(ind)=-isc(iabs(itype(i))) DU1(ind+1)=0.0d0 ind=ind+2 else DU1(ind)=mp/4-ip4 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 ! if (iabs(itype(i)).eq.ntyp1) cycle c write (iout,*) "i",i," itype",itype(i),ntyp1 if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) 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 enddo DMorig=DM DU1orig=DU1 DU2orig=DU2 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 #else dimen=(nct-nnt+1)+nside dimen1=(nct-nnt)+(nct-nnt+1) dimen3=dimen*3 write (iout,*) "Degrees_of_freedom",dimen3 #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 c write (iout,*) "dimen",dimen," dimen1",dimen1," dimen3",dimen3 c Zeroing out A and fricmat do i=1,dimen do j=1,dimen A(i,j)=0.0D0 enddo enddo c 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 c 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 c 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 c 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 (gmatout) 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,maxres2,maxres2,A) endif c 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 (gmatout) then write (iout,'(//a)') "Gmat" call matout(dimen,dimen,maxres2,maxres2,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 c Invert the G matrix call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob) if (gmatout) then write (iout,'(//a)') "Ginv" call matout(dimen,dimen,maxres2,maxres2,Ginv) endif #ifdef MPI if (nfgtasks.gt.1) then myginv_ng_count=maxres2*my_ng_count call MPI_Allgather(maxres2*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) c call MPI_Scatterv(ginv(1,1),nginv_counts(0), c & nginv_start(0),MPI_DOUBLE_PRECISION,ginv, c & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) c 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 c write (iout,*) "Master's chunk of ginv" c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,ginv) endif #endif if (osob) then write (iout,*) "The G matrix is singular." stop endif c 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(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec, & ierr,iwork) if (gmatout) then write (iout,'(//a)') & "Eigenvectors and eigenvalues of the G matrix" call eigout(dimen,dimen,maxres2,maxres2,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 #endif return end c------------------------------------------------------------------------------- SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B) implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision A(LM2,LM3),B(LM2) integer nc,nr,lm2,lm3,ka,kb,kc,n,i,j 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 c------------------------------------------------------------------------------- SUBROUTINE MATOUT(NC,NR,LM2,LM3,A) implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision A(LM2,LM3) integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j 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 c------------------------------------------------------------------------------- SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A) implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision A(LM2,LM3) integer nc,nr,lm2,lm3,n,ka,kb,kc,i,j 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 c------------------------------------------------------------------------------- SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A) implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision A(LM2,LM3) integer nc,nr,lm2,lm3,ka,kb,kc,i,j,n 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 c----------------------------------------------------------------------------- SUBROUTINE MATOUTR(N,A) c Prints the lower fragment of a symmetric matix implicit none integer n double precision a(n*(n+1)/2) integer i,j,k,nlim,jlim,jlim1 CHARACTER*6 LINE6 / '------' / CHARACTER*12 LINE12 / '------------' / double precision B(10) include 'COMMON.IOUNITS' DO 1 I=1,N,10 NLIM=MIN0(I+9,N) WRITE (IOUT,1000) (K,K=I,NLIM) WRITE (IOUT,1020) LINE6,(LINE12,K=I,NLIM) 1000 FORMAT (/7X,10(I5,2X)) 1020 FORMAT (A6,10A7) DO 2 J=I,N JLIM=MIN0(J,NLIM) JLIM1=JLIM-I+1 DO 3 K=I,JLIM 3 B(K-I+1)=A(J*(J-1)/2+K) WRITE (IOUT,1010) J,(B(K),K=1,JLIM1) 2 CONTINUE 1 CONTINUE 1010 FORMAT (I3,3X,10(F7.2)) RETURN END #ifdef FIVEDIAG c--------------------------------------------------------------------------- subroutine fivediagmult(n,DM,DU1,DU2,x,y) implicit none 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 c--------------------------------------------------------------------------- subroutine fivediaginv_mult(ndim,forces,d_a_vec) implicit none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.LAGRANGE.5diag' include 'COMMON.INTERACT' integer ndim double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim), & xsolv(ndim),d_a_vec(6*nres) integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev 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 rs(i)=forces(3*(i-1)+j) enddo call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv) ind=1 do i=innt,inct if (itype(i).eq.10)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 Conevert 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 if (itype(i).eq.10 .or. itype(i).eq.ntyp1) 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(:,2*nres)=0.0d0 if (nnt.gt.1) then accel(:,0)=accel(:,1) accel(:,1)=0.0d0 endif do ichain=2,nchain accel(:,chain_border1(1,ichain)-1)= & accel(:,chain_border1(1,ichain)) accel(:,chain_border1(1,ichain))=0.0d0 enddo do ichain=2,nchain 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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,dimen3) #endif return end #else c--------------------------------------------------------------------------- SUBROUTINE ginv_mult(z,d_a_tmp) implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierr,ierror #endif include 'COMMON.SETUP' include 'COMMON.TIME1' include 'COMMON.MD' include 'COMMON.LAGRANGE' double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00 &,time01,zcopy(dimen3) integer i,j,k,ind #ifdef MPI if (nfgtasks.gt.1) then if (fg_rank.eq.0) then c 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 c print *,"Processor",myrank," BROADCAST iorder in GINV_MULT" endif c write (2,*) "time00",time00 c write (2,*) "Before Scatterv" c call flush(2) c write (2,*) "Whole z (for FG master)" c do i=1,dimen c write (2,*) i,z(i) c enddo c call MPI_Barrier(FG_COMM,IERROR) time00=MPI_Wtime() call MPI_Scatterv(z,ng_counts(0),ng_start(0), & MPI_DOUBLE_PRECISION, & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) c write (2,*) "My chunk of z" do i=1,3*my_ng_count z(i)=zcopy(i) c write (2,*) i,z(i) enddo c write (2,*) "After SCATTERV" c call flush(2) c 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 c write (2,*) "time_scatter",time_scatter c write (2,*) "dimen",dimen," dimen3",dimen3," my_ng_count", c & my_ng_count c 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 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1, c & Ginv(i,j),z((j-1)*3+k+1), c & Ginv(i,j)*z((j-1)*3+k+1) c temp(ind)=temp(ind)+Ginv(i,j)*z((j-1)*3+k+1) temp(ind)=temp(ind)+Ginv(j,i)*z((j-1)*3+k+1) enddo enddo enddo time_ginvmult=time_ginvmult+MPI_Wtime()-time01 c write (2,*) "Before REDUCE" c call flush(2) c write (2,*) "z before reduce" c do i=1,dimen c write (2,*) i,temp(i) c 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 c write (2,*) "After REDUCE" c 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 c write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1 c call flush(2) c & Ginv(i,j),z((j-1)*3+k+1), c & 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) c d_a_tmp(ind)=d_a_tmp(ind) c & +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 c--------------------------------------------------------------------------- #ifdef GINV_MULT SUBROUTINE ginv_mult_test(z,d_a_tmp) implicit none include 'DIMENSIONS' include 'COMMON.LAGRANGE' double precision z(dimen),d_a_tmp(dimen) double precision ztmp(dimen/3),dtmp(dimen/3) c do i=1,dimen c d_a_tmp(i)=0.0d0 c do j=1,dimen c d_a_tmp(i)=d_a_tmp(i)+Ginv(i,j)*z(j) c enddo c enddo c c 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 #endif c--------------------------------------------------------------------------- SUBROUTINE fricmat_mult(z,d_a_tmp) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer IERROR #endif include 'COMMON.MD' include 'COMMON.LAGRANGE' include 'COMMON.IOUNITS' include 'COMMON.SETUP' include 'COMMON.TIME1' #ifndef LANG0 include 'COMMON.LANGEVIN' #else include 'COMMON.LANGEVIN.lang0' #endif double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00 &,time01,zcopy(dimen3) #ifdef MPI if (nfgtasks.gt.1) then if (fg_rank.eq.0) then c 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 c print *,"Processor",myrank," BROADCAST iorder in FRICMAT_MULT" endif c call MPI_Barrier(FG_COMM,IERROR) time00=MPI_Wtime() call MPI_Scatterv(z,ng_counts(0),ng_start(0), & MPI_DOUBLE_PRECISION, & zcopy,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) c write (2,*) "My chunk of z" do i=1,3*my_ng_count z(i)=zcopy(i) c 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)*z((j-1)*3+k+1) enddo enddo enddo time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01 c write (2,*) "Before REDUCE" c write (2,*) "d_a_tmp before reduce" c do i=1,dimen3 c write (2,*) i,temp(i) c enddo c 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 c write (2,*) "After REDUCE" c 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 c write (iout,*) "Vector d_a" c do i=1,dimen3 c write (2,*) i,d_a_tmp(i) c enddo return end #endif