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)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
+ integer time00
#endif
include 'COMMON.VAR'
include 'COMMON.CHAIN'
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'
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
& 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
+ end
c------------------------------------------------------------------
subroutine setup_MD_matrices
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
- integer ierror
+ integer ierror,ierr
+ double precision time00
#endif
+ include 'COMMON.CONTROL'
include 'COMMON.SETUP'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
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
+ 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),Gcopy(maxres2,maxres2),
- & Ghalf(mmaxres2),sqreig(maxres2)
+ 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()
A(k+ii,jj)=1.0d0
enddo
enddo
- if (lprn) then
+ if (gmatout) then
write (iout,*)
write (iout,*) "Vector massvec"
do i=1,dimen1
enddo
enddo
- if (lprn) then
+ if (gmatout) then
write (iout,'(//a)') "Gmat"
call matout(dimen,dimen,maxres2,maxres2,Gmat)
endif
enddo
c Invert the G matrix
call MATINVERT(dimen,maxres2,Gcopy,Ginv,osob)
- if (lprn) then
+ if (gmatout) then
write (iout,'(//a)') "Ginv"
call matout(dimen,dimen,maxres2,maxres2,Ginv)
endif
enddo
call gldiag(maxres2,dimen,dimen,Ghalf,work,Geigen,Gvec,
& ierr,iwork)
- if (lprn) then
+ if (gmatout) then
write (iout,'(//a)')
& "Eigenvectors and eigenvalues of the G matrix"
call eigout(dimen,dimen,maxres2,maxres2,Gvec,Geigen)
enddo
enddo
endif
+#endif
return
end
c-------------------------------------------------------------------------------
SUBROUTINE EIGOUT(NC,NR,LM2,LM3,A,B)
- implicit real*8 (a-h,o-z)
+ 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)
END
c-------------------------------------------------------------------------------
SUBROUTINE MATOUT(NC,NR,LM2,LM3,A)
- implicit real*8 (a-h,o-z)
+ 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)
END
c-------------------------------------------------------------------------------
SUBROUTINE MATOUT1(NC,NR,LM2,LM3,A)
- implicit real*8 (a-h,o-z)
+ 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)
END
c-------------------------------------------------------------------------------
SUBROUTINE MATOUT2(NC,NR,LM2,LM3,A)
- implicit real*8 (a-h,o-z)
+ 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)
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 real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
- integer ierr
+ 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---------------------------------------------------------------------------
#ifdef GINV_MULT
SUBROUTINE ginv_mult_test(z,d_a_tmp)
+ implicit none
include 'DIMENSIONS'
- integer dimen
-c include 'COMMON.MD'
+ include 'COMMON.LAGRANGE'
double precision z(dimen),d_a_tmp(dimen)
double precision ztmp(dimen/3),dtmp(dimen/3)
integer IERROR
#endif
include 'COMMON.MD'
+ include 'COMMON.LAGRANGE'
include 'COMMON.IOUNITS'
include 'COMMON.SETUP'
include 'COMMON.TIME1'
c enddo
return
end
+#endif
projects / unres.git / blobdiff