use comm_cipiszcze
use energy_data
use geometry_data, only: nres
- use control_data
+ use control_data !el, only: mucadyn,lmuca
#ifdef MPI
include 'mpif.h'
real(kind=8) :: time00
! include 'COMMON.CONTROL'
! include 'COMMON.MUCA'
! include 'COMMON.TIME1'
-
- integer :: i,j,ind,itime
+ integer ::i,j,ind,itime,k
real(kind=8) :: zapas(6*nres) !,muca_factor !maxres6=6*maxres
+ real(kind=8) :: rs(dimen),xsolv(dimen)
+#ifdef CHECK5DSOL
+ real(kind=8) :: rscheck(dimen),rsold(dimen)
+#endif
logical :: lprn = .false.
!el common /cipiszcze/ itime
itime = itt_comm
#ifdef TIMING
time00=MPI_Wtime()
#endif
+#ifdef FIVEDIAG
+ if (lprn) then
+ write (iout,*) "Potential forces backbone"
+ do i=nnt,nct
+ 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 j=1,3
+ ind=1
+ do i=nnt,nct
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1)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"
+ do i=1,dimen
+ write(iout,*) i,rs(i)
+ enddo
+ endif
+
+ call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
+ if (lprn) then
+ write (iout,*) "Solution of the 5-diagonal equations system"
+ do i=1,dimen
+ write (iout,'(i5,f10.5)') i,xsolv(i)
+ enddo
+ endif
+#ifdef CHECK5DSOL
+! Check the solution
+ rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
+ DU2orig(1)*xsolv(3)
+ rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
+ DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
+
+ do i=3,dimen-2
+ rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
+ xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
+ xsolv(i+1)+DU2orig(i)*xsolv(i+2)
+ enddo
+ rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
+ DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
+ xsolv(dimen-1)+DU1orig(dimen-1)*&
+ xsolv(dimen)
+ rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
+ xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
+
+ do i=1,dimen
+ write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
+ "ratio",rscheck(i)/rsold(i)
+ enddo
+! end check
+#endif
+ ind=1
+ do i=nnt,nct
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1)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
+ if (lprn) then
+ do i=nnt,nct
+ 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
+ do j=1,3
+ d_a(j,0)=d_a(j,nnt)
+ enddo
+ do i=nnt,nct
+ 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
+
+ if (lprn) then
+ write(iout,*) 'acceleration 3D FIVEDIAG'
+ 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
+#else
+! Old procedure
do j=1,3
zapas(j)=-gcart(j,0)
enddo
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)
+ i,(-gxcart(j,i),j=1,3)
do j=1,3
ind=ind+1
zapas(ind)=-gxcart(j,i)
enddo
endif
enddo
-
+
+#endif
if(lmuca) then
imtime=imtime+1
if(mucadyn.gt.0) call muca_update(potE)
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.TIME1'
logical :: lprn = .false.
logical :: osob
+#ifndef FIVEDIAG
+ real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
+#endif
real(kind=8) :: dtdi
real(kind=8),dimension(2*nres) :: 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)
integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres
!el common /przechowalnia/ Gcopy,Ghalf
real(kind=8) :: coeff
- integer :: i,j,ind,ind1,k,ii,jj,m,m1,ii1,iti,nres2,ierr
+ integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
+ real(kind=8) :: ip4
+ integer :: iz
+ relfeh=1.0d-14
nres2=2*nres
+ write (iout,*) "before FIVEDIAG"
+#ifndef FIVEDIAG
+ write (iout,*) "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)
+#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)
dimen=(nct-nnt+1)+nside
dimen1=(nct-nnt)+(nct-nnt+1)
dimen3=dimen*3
+ write (iout,*) "nnt",nnt," nct",nct," nside",nside
+#ifdef FIVEDIAG
+ 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
+! if (iabs(itype(i)).eq.ntyp1) cycle
+ DM(ind)=2*ip4+mp/2
+ if (iabs(itype(i)).eq.10 .or. iabs(itype(i)).eq.ntyp1) then
+ if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10)
+ ind=ind+1
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(i)))
+ DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
+ 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
+ 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
+
+ ind=1
+ do i=nnt,nct-1
+! if (iabs(itype(i)).eq.ntyp1) cycle
+ 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
+ DMorig=DM
+ DU1orig=DU1
+ DU2orig=DU2
+ write (iout,*) "nind",nind," dimen",dimen
+ if (nind.ne.dimen) then
+ write (iout,*) "ERROR, dimensions differ"
+#ifdef MPI
+ call MPI_Finalize(ierr)
+#endif
+ stop
+ endif
+ write (iout,*) "The upper part of the five-diagonal inertia matrix"
+ do i=1,dimen
+ if (i.lt.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
+ else if (i.eq.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+ else
+ write (iout,'(i3,3f10.5)') i,DM(i)
+ endif
+ enddo
+
+ call FDISYP (dimen, DM, DU1, DU2, MARK)
+
+ if (mark.eq.-1) then
+ write (iout,*) "ERROR: the inertia matrix is not positive definite"
+#ifdef MPI
+ call MPI_Finalize(ierr)
+#endif
+ stop
+ else if (mark.eq.0) then
+ write (iout,*) "ERROR: the inertia matrix is singular"
+#ifdef MPI
+ call MPI_Finalize(ierr)
+#endif
+ else if (mark.eq.1) then
+ write (iout,*) "The transformed inertia matrix"
+ do i=1,dimen
+ if (i.lt.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
+ else if (i.eq.dimen-1) then
+ write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+ else
+ write (iout,'(i3,3f10.5)') i,DM(i)
+ endif
+ enddo
+ endif
+! Diagonalization of the pentadiagonal matrix
+ allocate(DDU1(2*nres))
+ allocate(DDU2(2*nres))
+ allocate(DDM(2*nres))
+ do i=1,dimen-1
+ DDU1(i+1)=DU1orig(i)
+ enddo
+ do i=1,dimen-2
+ DDU2(i+2)=DU2orig(i)
+ enddo
+ DDM=DMorig
+ call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
+ if (lprn) then
+ write (iout,*) &
+ "Eigenvalues of the five-diagonal inertia matrix"
+ do i=1,dimen
+ write (iout,'(i5,f10.5)') i,geigen(i)
+ enddo
+ endif
+ deallocate(DDU1)
+ deallocate(DDU2)
+ deallocate(DDM)
+#else
+! Old Gmatrix
#ifdef MPI
if (nfgtasks.gt.1) then
time00=MPI_Wtime()
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
+ if (itype(i).eq.10) 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
+ do i=1,dimen
+ do j=1,dimen
+ do k=1,dimen
+ do l=1,dimen
+ Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
+ enddo
+ enddo
+ enddo
+ enddo
+ 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
enddo
endif
deallocate(Gcopy)
+#endif
return
end subroutine setup_MD_matrices
!-----------------------------------------------------------------------------
projects / unres4.git / blobdiff