! include 'COMMON.CONTROL'
! include 'COMMON.MUCA'
! include 'COMMON.TIME1'
- integer ::i,j,ind,itime,k
+ integer ::i,j,ind,itime,mnum
real(kind=8) :: zapas(6*nres) !,muca_factor !maxres6=6*maxres
real(kind=8) :: rs(dimen),xsolv(dimen)
#ifdef CHECK5DSOL
enddo
write (iout,*) "Potential forces sidechain"
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.eq.5)&
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
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.eq.5)then
rs(ind)=-gcart(j,i)-gxcart(j,i)
ind=ind+1
else
#endif
ind=1
do i=nnt,nct
- if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum))then
d_a(j,i)=xsolv(ind)
ind=ind+1
else
d_a(j,0)=d_a(j,nnt)
enddo
do i=nnt,nct
- if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+ mnum=molnum(i)
+! if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum))then
do j=1,3
d_a(j,i)=d_a(j,i+1)-d_a(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
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.ne.5) then
if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
i,(-gxcart(j,i),j=1,3)
do j=1,3
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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.ne.5) then
do j=1,3
ind=ind+1
d_a(j,i+nres)=d_a_work(ind)
logical :: lprn = .false.
logical :: osob
#ifndef FIVEDIAG
- real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
+ real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult
#endif
real(kind=8) :: dtdi
real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres
real(kind=8),dimension(8*6*nres) :: work !(8*maxres6)
integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres
!el common /przechowalnia/ Gcopy,Ghalf
- real(kind=8) :: coeff
+ 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
+ integer :: iz,mnum
+ print *,"just entered"
relfeh=1.0d-14
nres2=2*nres
-
+ print *,"FIVEDIAG"
write (iout,*) "before FIVEDIAG"
#ifndef FIVEDIAG
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
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
nind=2
endif
do i=nnt+1,nct-1
-! if (iabs(itype(i)).eq.ntyp1) cycle
+! if (iabs(itype(i,1)).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)
+ if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) 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)))
- DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
+ 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
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)))
+ mnum=molnum(i)
+ if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
+ .and.mnum.eq.5) then
+ DU1(ind)=-isc(iabs(itype(i,1)))
DU1(ind+1)=0.0d0
ind=ind+2
else
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
+ 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.ne.5) then
DU2(ind)=mp/4-ip4
DU2(ind+1)=0.0d0
ind=ind+2
ind=ind+1
endif
enddo
+#else
+ ip4=ip(1)/4
+ DM(1)=mp(1)/4+ip4
+ if (iabs(itype(nnt,1)).eq.10) then
+ DM(1)=DM(1)+msc(10,1)
+ ind=2
+ nind=1
+ else
+ DM(1)=DM(1)+isc(iabs(itype(nnt,1)),1)
+ DM(2)=msc(iabs(itype(nnt,1)),1)+isc(iabs(itype(nnt,1)),1)
+ ind=3
+ nind=2
+ endif
+ do i=nnt+1,nct-1
+! if (iabs(itype(i,1)).eq.ntyp1) cycle
+ DM(ind)=2*ip4+mp(1)/2
+ if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
+ if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10,1)
+ ind=ind+1
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(i,1)),1)
+ DM(ind+1)=msc(iabs(itype(i,1)),1)+isc(iabs(itype(i,1)),1)
+ ind=ind+2
+ endif
+ enddo
+ if (nct.gt.nnt) then
+ DM(ind)=ip4+mp(1)/4
+ if (iabs(itype(nct,1)).eq.10) then
+ DM(ind)=DM(ind)+msc(10,1)
+ nind=ind
+ else
+ DM(ind)=DM(ind)+isc(iabs(itype(nct,1)),1)
+ DM(ind+1)=msc(iabs(itype(nct,1)),1)+isc(iabs(itype(nct,1)),1)
+ 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,mnum)).ne.ntyp1 &
+ .and.mnum.eq.5) then
+ DU1(ind)=-isc(iabs(itype(i,1)),1)
+ DU1(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU1(ind)=mp(1)/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.ne.5) then
+ DU2(ind)=mp(1)/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
+#endif
DMorig=DM
DU1orig=DU1
DU2orig=DU2
! 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
- massvec(ind1)=mp
+ coeff=0.25d0*IP(mnum)
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ print *,i,coeff,mp(mnum)
+ massvec(ind1)=mp(mnum)
Gmat(ind,ind)=coeff
+ print *,"i",mp(mnum)
A(ind1,ind)=0.5d0
enddo
m1=nct-nnt+1
ind=0
ind1=0
- msc(ntyp1)=1.0d0
+ do i=1,2
+ msc(ntyp1_molec(i),i)=1.0d0
+ enddo
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
+ 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.ne.5) then
ind1=ind1+1
ii1= ind1+m1
A(ii,ii1)=1.0d0
- Gmat(ii1,ii1)=ISC(iabs(iti))
+ 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
- iti=itype(i)
+ iti=itype(i,1)
ind=ind+1
do j=nnt,i
ii = ind
Bmat(1,1)=1.0d0
j=2
do i=nnt,nct
- if (itype(i).eq.10) then
+ mnum=molnum(i)
+ if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
+ .or.(mnum.eq.5)) then
Bmat(i-nnt+2,j-1)=-1
Bmat(i-nnt+2,j)=1
j=j+1
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)
projects / unres4.git / blobdiff