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 real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #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 double precision zapas(MAXRES6),muca_factor logical lprn /.false./ common /cipiszcze/ itime #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) 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) 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 #ifdef TIMING time_lagrangian=time_lagrangian+MPI_Wtime()-time00 #endif return end c------------------------------------------------------------------ subroutine setup_MD_matrices implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierror #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' 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' integer i,j logical lprn /.false./ logical osob double precision dtdi,massvec(maxres2),Gcopy(maxres2,maxres2), & Ghalf(mmaxres2),sqreig(maxres2), invsqreig(maxres2), Gcopytmp, & Gsqrptmp, Gsqrmtmp, Gvec2tmp,Gvectmp(maxres2,maxres2) double precision work(8*maxres6) integer iwork(maxres6) common /przechowalnia/ Gcopy,Ghalf,invsqreig,Gvectmp 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) 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 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 do i=nnt,nct ind=ind+1 ii = ind+m iti=itype(i) massvec(ii)=msc(iti) if (iti.ne.10) then ind1=ind1+1 ii1= ind1+m1 A(ii,ii1)=1.0d0 Gmat(ii1,ii1)=ISC(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 (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,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 (lprn) 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 (lprn) 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) if (lprn .and. (me.eq.king .or. .not. out1file) ) then 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) endif 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 (lprn) 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)) invsqreig(i)=1.d0/sqreig(i) enddo do i=1,dimen do j=1,dimen Gvectmp(i,j)=Gvec(j,i) enddo enddo do i=1,dimen do j=1,dimen Gsqrptmp=0.0d0 Gsqrmtmp=0.0d0 Gcopytmp=0.0d0 do k=1,dimen c Gvec2tmp=Gvec(i,k)*Gvec(j,k) Gvec2tmp=Gvec(k,i)*Gvec(k,j) Gsqrptmp=Gsqrptmp+Gvec2tmp*sqreig(k) Gsqrmtmp=Gsqrmtmp+Gvec2tmp*invsqreig(k) Gcopytmp=Gcopytmp+Gvec2tmp*Geigen(k) enddo Gsqrp(i,j)=Gsqrptmp Gsqrm(i,j)=Gsqrmtmp Gcopy(i,j)=Gcopytmp enddo enddo do i=1,dimen do j=1,dimen Gvec(i,j)=Gvectmp(j,i) 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 return end c------------------------------------------------------------------------------- SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision 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 c------------------------------------------------------------------------------- SUBROUTINE MATOUT(NC,NR,LM2,LM3,A) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision 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 c------------------------------------------------------------------------------- SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision 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 c------------------------------------------------------------------------------- SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' double precision 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 c--------------------------------------------------------------------------- SUBROUTINE ginv_mult(z,d_a_tmp) implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierr #endif include 'COMMON.SETUP' include 'COMMON.TIME1' include 'COMMON.MD' double precision z(dimen3),d_a_tmp(dimen3),temp(maxres6),time00 &time01 #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, & z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) c write (2,*) "My chunk of z" c do i=1,3*my_ng_count c write (2,*) i,z(i) c 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) include 'DIMENSIONS' integer dimen c include 'COMMON.MD' 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.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 #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, & z,3*my_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR) c write (2,*) "My chunk of z" c do i=1,3*my_ng_count c write (2,*) i,z(i) c 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