! Groups the ALPHAs and the BETAs
do k=1,numch
do j=2,nres-1
- if(itype(j).ne.10) then
+ if(itype(j,1).ne.10) then
indg=indg+1
ngroup(indg)=2
do i=1,ngroup(indg)
#endif
endif
do j=jstart,jend
- if (itype(j).eq.10) then
+ if (itype(j,1).eq.10) then
iang=2
else
iang=4
enddo
if (ishift.gt.0) then
do j=0,ishift-1
- if (itype(jend+j).eq.10) then
+ if (itype(jend+j,1).eq.10) then
iang=2
else
iang=4
enddo
else
do j=0,-ishift-1
- if (itype(jstart+j).eq.10) then
+ if (itype(jstart+j,1).eq.10) then
iang=2
else
iang=4
! enddo
! call intout
do j=jstart,jend
- if (itype(j).eq.10) then
+ if (itype(j,1).eq.10) then
iang=2
else
iang=4
do i=1,nwindow
i1=winstart(i)
i2=winend(i)
- it1=itype(i1)
- it2=itype(i2)
- write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,&
+ it1=itype(i1,1)
+ it2=itype(i2,1)
+ write (iout,'(a,i3,a,i3,a,i3)') restyp(it1,1),i1,restyp(it2,1),i2,&
' length',winlen(i)
enddo
endif
! print *,'nnt=',nnt,' nct=',nct
ngly=0
do i=nnt,nct
- if (itype(i).eq.10) ngly=ngly+1
+ if (itype(i,1).eq.10) ngly=ngly+1
enddo
mmm=nct-nnt-ngly+1
if (mmm.gt.0) then
error=.true.
return
endif
- if (itype(inds).eq.10) goto 111
+ if (itype(inds,1).eq.10) goto 111
do j=1,i-1
if (inds.eq.ind_side(j)) goto 111
enddo
! Carry out perturbation.
do i=1,nside_move
ii=ind_side(i)
- iti=itype(ii)
+ iti=itype(ii,1)
call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
if (fail) then
isctry=isctry+1
goto 301
endif
if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') &
- 'Side chain ',restyp(iti),ii,' moved to ',&
+ 'Side chain ',restyp(iti,1),ii,' moved to ',&
alph(ii)*rad2deg,omeg(ii)*rad2deg
enddo
moves(3)=moves(3)+1
!------------------------------------------------------------------------------
! THETA move
400 end_select=iran_num(3,nres)
- theta_new=gen_theta(itype(end_select),phi(end_select),&
+ theta_new=gen_theta(itype(end_select,1),phi(end_select),&
phi(end_select+1))
if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') &
'Theta ',end_select,' moved from',theta(end_select)*rad2deg,&
nbond=nct-nnt
do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
+ if (itype(i,1).ne.10) nbond=nbond+1
enddo
!
if (lprn1) then
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind1=ind1+1
do j=1,3
Bmat(ind+j,ind1)=dC_norm(j,i+nres)
Td(i)=Td(i)+vbl*Tmat(i,ind)
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
- Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
+ Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
endif
enddo
enddo
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
do j=1,3
ind=ind+1
zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
i,(dC(j,i),j=1,3),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),xx
endif
endif
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
- ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
- diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
+ ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
+ diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
if (diffbond.gt.diffmax) diffmax=diffbond
if (ppvec(ind).gt.0.0d0) then
ppvec(ind)=dsqrt(ppvec(ind))
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
do j=1,3
dc(j,i+nres)=zapas(ind+j)
dc_work(ind+j)=zapas(ind+j)
i,(dC(j,i),j=1,3),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),xx
endif
#endif
KEt_p=0.0d0
KEt_sc=0.0d0
-! write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
+! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
! The translational part for peptide virtual bonds
do j=1,3
incr(j)=d_t(j,0)
incr(j)=d_t(j,0)
enddo
do i=nnt,nct
- iti=iabs(itype(i))
- if (itype(i).eq.10) then
+ iti=iabs(itype(i,1))
+ if (itype(i,1).eq.10) then
do j=1,3
v(j)=incr(j)
enddo
! The rotational part of the side chain virtual bond
KEr_sc=0.0D0
do i=nnt,nct
- iti=iabs(itype(i))
- if (itype(i).ne.10) then
+ iti=iabs(itype(i,1))
+ if (itype(i,1).ne.10) then
do j=1,3
incr(j)=d_t(j,nres+i)
enddo
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
ind=ind+1
v_work(ind)=d_t(j,i+nres)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
adt=d_a_old(j,inres)*d_time
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
enddo
endif
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
enddo
endif
! Side chains
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
epdriftij= &
dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=fact*d_t(j,inres)
stdforcp(i)=stdfp*dsqrt(gamp)
enddo
do i=nnt,nct
- stdforcsc(i)=stdfsc(iabs(itype(i))) &
- *dsqrt(gamsc(iabs(itype(i))))
+ stdforcsc(i)=stdfsc(iabs(itype(i,1))) &
+ *dsqrt(gamsc(iabs(itype(i,1))))
enddo
endif
! Open the pdb file for snapshotshots
Ek1=0.0d0
ii=0
do i=nnt,nct
- if (itype(i).eq.10) then
+ if (itype(i,1).eq.10) then
j=ii+3
else
j=ii+6
0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
enddo
endif
- if (itype(i).ne.10) ii=ii+3
- write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i))
+ if (itype(i,1).ne.10) ii=ii+3
+ write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1))
write (iout,*) "ii",ii
do k=1,3
ii=ii+1
write (iout,*) "k",k," ii",ii,"EK1",EK1
- if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2
- Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2
+ if (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1)))*(d_t_work(ii)-d_t_work(ii-3))**2
+ Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)))*d_t_work(ii)**2
enddo
write (iout,*) "i",i," ii",ii
enddo
d_t(j,0)=d_t(j,nnt)
enddo
do i=nnt,nct
- if (itype(i).eq.10) then
+ if (itype(i,1).eq.10) then
do j=1,3
d_t(j,i)=d_t(j,i+1)-d_t(j,i)
enddo
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
ind=ind+1
d_t(j,i+nres)=d_t_work(ind)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
dc_work(ind+j)=dc_old(j,i+nres)
d_t_work(ind+j)=d_t_old(j,i+nres)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
inres=i+nres
do j=1,3
dc(j,inres)=dc_work(ind+j)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_work(ind+j)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
dc_work(ind+j)=dc_old(j,i+nres)
d_t_work(ind+j)=d_t_old(j,i+nres)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
dc(j,inres)=dc_work(ind+j)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_work(ind+j)
enddo
M_SC=0.0d0
do i=nnt,nct
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
M_SC=M_SC+msc(iabs(iti))
inres=i+nres
do j=1,3
enddo
do i=nnt,nct
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
inres=i+nres
do j=1,3
pr(j)=c(j,inres)-cm(j)
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
- iti=iabs(itype(i))
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+ iti=iabs(itype(i,1))
inres=i+nres
Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
dc_norm(1,inres))*vbld(inres)*vbld(inres)
enddo
enddo
do i=nnt,nct
- if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
call vecpr(vrot(1),dc(1,inres),vp)
do j=1,3
incr(j)=d_t(j,0)
enddo
do i=nnt,nct
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
inres=i+nres
do j=1,3
pr(j)=c(j,inres)-cm(j)
enddo
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
v(j)=incr(j)+d_t(j,inres)
enddo
L(j)=L(j)+msc(iabs(iti))*vp(j)
enddo
! write (iout,*) "L",(l(j),j=1,3)
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
v(j)=incr(j)+d_t(j,inres)
enddo
vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
enddo
endif
- amas=msc(iabs(itype(i)))
+ amas=msc(iabs(itype(i,1)))
summas=summas+amas
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
enddo
if (lprn) write (iout,*) "RATTLE1"
nbond=nct-nnt
do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
+ if (itype(i,1).ne.10) nbond=nbond+1
enddo
! Make a folded form of the Ginv-matrix
ind=0
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ii=ii+1
do j=1,3
ind=ind+1
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind1=ind1+1
do j=1,3
dC_uncor(j,ind1)=dC(j,i+nres)
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
x(ind)=vbld(i+1)**2-vbl**2
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
endif
enddo
if (lprn) then
i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
do j=1,3
xx=0.0d0
i,(dC(j,i),j=1,3),x(ind),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),x(ind),xx
endif
i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
do j=1,nbond
Cmat(ind,j)=0.0d0
x(ind)=scalar(d_t(1,i),dC(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
endif
i,ind,(d_t(j,i),j=1,3),x(ind)
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
do j=1,3
xx=0.0d0
i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,2e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
if (lprn) write (iout,*) "RATTLE_BROWN"
nbond=nct-nnt
do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
+ if (itype(i,1).ne.10) nbond=nbond+1
enddo
! Make a folded form of the Ginv-matrix
ind=0
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ii=ii+1
do j=1,3
ind=ind+1
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind1=ind1+1
do j=1,3
dC_uncor(j,ind1)=dC(j,i+nres)
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
x(ind)=vbld(i+1)**2-vbl**2
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
endif
enddo
if (lprn) then
i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),&
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
do j=1,3
xx=0.0d0
i,(dC(j,i),j=1,3),x(ind),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),x(ind),xx
endif
i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
ind=ind+3
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
d_t_work(ind+j)=d_t(j,i+nres)
enddo
ind=ind+3
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
friction(j,i+nres)=fric_work(ind+j)
enddo
do i=nnt,nct
ind=ind+1
ii = ind+m
- iti=itype(i)
+ iti=itype(i,1)
gamvec(ii)=gamsc(iabs(iti))
enddo
if (surfarea) call sdarea(gamvec)
do j=1,3
ff(j)=ff(j)+force(j,i)
enddo
- if (itype(i+1).ne.ntyp1) then
+ if (itype(i+1,1).ne.ntyp1) then
do j=1,3
stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
ff(j)=ff(j)+force(j,i+nres+1)
stochforc(j,0)=ff(j)+force(j,nnt+nres)
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
stochforc(j,i+nres)=force(j,i+nres)
enddo
ind=ind+3
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
stochforcvec(ind+j)=stochforc(j,i+nres)
enddo
enddo
! Load side chain radii
do i=nnt,nct
- iti=itype(i)
+ iti=itype(i,1)
radius(i+nres)=restok(iti)
enddo
! do i=1,2*nres
ind=ind+1
gamvec(ind) = ratio * gamvec(ind)
enddo
- stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
+ stdforcsc(i)=stdfsc(itype(i,1))*dsqrt(gamvec(ind))
if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
endif
enddo
enddo
write (iout,*) "Potential forces sidechain"
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+ if (itype(i,1).ne.10 .and. itype(i,1).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
+ if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)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
+ if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)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
+ if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) 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
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) 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
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
ind=ind+1
d_a(j,i+nres)=d_a_work(ind)
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)))
+ if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) 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
+! 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,1)).ne.ntyp1) then
DU2(ind)=mp/4-ip4
DU2(ind+1)=0.0d0
ind=ind+2
do i=nnt,nct
ind=ind+1
ii = ind+m
- iti=itype(i)
+ iti=itype(i,1)
massvec(ii)=msc(iabs(iti))
if (iti.ne.10 .and. iti.ne.ntyp1) then
ind1=ind1+1
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
+ if (itype(i,1).eq.10) then
Bmat(i-nnt+2,j-1)=-1
Bmat(i-nnt+2,j)=1
j=j+1
!el local variables
integer :: it,i
- it=itype(2)
+ it=itype(2,1)
do i=1,101
vbld(nres+2)=0.5d0+0.05d0*(i-1)
call chainbuild
ncont=0
kkk=3
do i=nnt+kkk,nct
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
do j=nnt,i-kkk
- itj=iabs(itype(j))
+ itj=iabs(itype(j,1))
if (ipot.ne.4) then
! rcomp=sigmaii(iti,itj)+1.0D0
rcomp=facont*sigmaii(iti,itj)
do i=1,ncont
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,1)
+ it2=itype(i2,1)
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,1),i1,restyp(it2,1),i2
enddo
endif
co = 0.0d0
do i=1,ncont
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,1)
+ it2=itype(i2,1)
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,1),i1,restyp(it2,1),i2
enddo
endif
! finding hairpins
ii1=iharp(3,i)
jj1=iharp(4,i)
write (iout,*)
- write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1)
- write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
+ write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1)
+ write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1)
! do k=jj1,j1,-1
-! write (iout,'(a,i3,$)') restyp(itype(k)),k
+! write (iout,'(a,i3,$)') restyp(itype(k,1)),k
! enddo
enddo
endif
ees=0.0
evdw=0.0
do 1 i=nnt,nct-2
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
xi=c(1,i)
yi=c(2,i)
zi=c(3,i)
ymedi=yi+0.5*dyi
zmedi=zi+0.5*dzi
do 4 j=i+2,nct-1
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
iteli=itel(i)
itelj=itel(j)
if (j.eq.i+2 .and. itelj.eq.2) iteli=2
do i=1,ncont
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,1)
+ it2=itype(i2,1)
write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
- i,restyp(it1),i1,restyp(it2),i2,econt(i)
+ i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
enddo
endif
! For given residues keep only the contacts with the greatest energy.
do i=1,ncont
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,1)
+ it2=itype(i2,1)
write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
- i,restyp(it1),i1,restyp(it2),i2,econt(i)
+ i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
enddo
endif
return
do kkk=1,nperm
! do i=nz_start,nz_end
! iatom=iatom+1
-! iti=itype(i)
+! iti=itype(i,1)
! do k=1,3
! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
! crefcopy(k,iatom,kkk)=cref(k,i,kkk)
iatom=0
do i=nz_start,nz_end
iatom=iatom+1
- iti=itype(i)
+ iti=itype(i,1)
do k=1,3
ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
crefcopy(k,iatom)=cref(k,i,kkk)
iatom=0
do i=nz_start,nz_end
iatom=iatom+1
- iti=itype(i)
+ iti=itype(i,1)
do k=1,3
ccopy(k,iatom)=c(k,i)
crefcopy(k,iatom)=crefjlee(k,i)
c(j,i)=cart_base(j,i+ist,ii)
! cref(j,i)=c(j,i)
enddo
-!d write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3)
+!d write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3)
enddo
!d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
!d non_conv)
!d write (iout,'(a,f10.5)')
!d & 'Initial RMS deviation from reference structure:',rms
- if (itype(nres).eq.ntyp1) then
+ if (itype(nres,1).eq.ntyp1) then
do j=1,3
dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
endif
- if (itype(1).eq.ntyp1) then
+ if (itype(1,1).eq.ntyp1) then
do j=1,3
dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
link_end0=link_end
link_end=min0(link_end,nss)
do i=nnt,nct
- if (itype(i).ne.10) then
-!d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
- call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
+ if (itype(i,1).ne.10) then
+!d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
+ call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
endif
enddo
call chainbuild
glycine=.true.
do while(glycine)
ind_sc=iran_num(nnt,nct)
- glycine=(itype(ind_sc).eq.10)
+ glycine=(itype(ind_sc,1).eq.10)
enddo
alph0=alph(ind_sc)
omeg0=omeg(ind_sc)
- call gen_side(itype(ind_sc),theta(ind_sc+1),alph(ind_sc),&
+ call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
omeg(ind_sc),fail)
call chainbuild
call etotal(energia)
! enddo !iblock
! do i=1,maxres
-! itype(i)=0
+! itype(i,1)=0
! itel(i)=0
! enddo
! Initialize the bridge arrays
ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,&
ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,&
ichunk,int_index_old
-
+ integer,dimension(5) :: nct_molec
!el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
!el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
ii=nint_gr(i)+1
call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,&
- iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
+ iatsc_s,iatsc_e,jj+1,nct_molec(1),nint_gr(i),istart(i,ii),iend(i,ii),*12)
#else
nint_gr(i)=2
istart(i,1)=i+1
iend(i,1)=jj-1
istart(i,2)=jj+1
- iend(i,2)=nct
+ iend(i,2)=nct_molec(1)
#endif
endif
else
#ifdef MPI
call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,&
- iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
+ iatsc_s,iatsc_e,i+1,nct_molec(1),nint_gr(i),istart(i,1),iend(i,1),*12)
#else
nint_gr(i)=1
istart(i,1)=i+1
- iend(i,1)=nct
- ind_scint=ind_scint+nct-i
+ iend(i,1)=nct_molec(1)
+ ind_scint=ind_scint+nct_molec(1)-i
#endif
endif
#ifdef MPI
ispp=4 !?? wham ispp=2
#ifdef MPI
! Now partition the electrostatic-interaction array
- npept=nct-nnt
+ if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+ npept=nres_molec(1)-nnt-1
+ else
+ npept=nres_molec(1)-nnt
+ endif
nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
if (lprint) &
iatel_e=0
ind_eleint=0
ind_eleint_old=0
- do i=nnt,nct-3
+ if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+ nct_molec(1)=nres_molec(1)-1
+ else
+ nct_molec(1)=nres_molec(1)
+ endif
+! print *,"nct",nct,nct_molec(1),itype(nres_molec(1),1),ntyp_molec(1)
+ do i=nnt,nct_molec(1)-3
ijunk=0
call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,&
- iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
+ iatel_s,iatel_e,i+ispp,nct_molec(1)-1,ijunk,ielstart(i),ielend(i),*13)
enddo ! i
13 continue
if (iatel_s.eq.0) iatel_s=1
ind_eleint_vdw_old=0
iatel_s_vdw=0
iatel_e_vdw=0
- do i=nnt,nct-3
+ do i=nnt,nct_molec(1)-3
ijunk=0
call int_partition(ind_eleint_vdw,my_ele_inds_vdw,&
my_ele_inde_vdw,i,&
- iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),&
+ iatel_s_vdw,iatel_e_vdw,i+2,nct_molec(1)-1,ijunk,ielstart_vdw(i),&
ielend_vdw(i),*15)
! write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
! & " ielend_vdw",ielend_vdw(i)
15 continue
#else
iatel_s=nnt
- iatel_e=nct-5 ! ?? wham iatel_e=nct-3
+ iatel_e=nct_molec(1)-5 ! ?? wham iatel_e=nct-3
do i=iatel_s,iatel_e
ielstart(i)=i+4 ! ?? wham +2
- ielend(i)=nct-1
+ ielend(i)=nct_molec(1)-1
enddo
iatel_s_vdw=nnt
- iatel_e_vdw=nct-3
+ iatel_e_vdw=nct_molec(1)-3
do i=iatel_s_vdw,iatel_e_vdw
ielstart_vdw(i)=i+2
- ielend_vdw(i)=nct-1
+ ielend_vdw(i)=nct_molec(1)-1
enddo
#endif
if (lprint) then
iatscp_e=0
ind_scpint=0
ind_scpint_old=0
- do i=nnt,nct-1
+ do i=nnt,nct_molec(1)-1
if (i.lt.nnt+iscp) then
!d write (iout,*) 'i.le.nnt+iscp'
call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,&
- iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),&
+ iatscp_s,iatscp_e,i+iscp,nct_molec(1),nscp_gr(i),iscpstart(i,1),&
iscpend(i,1),*14)
else if (i.gt.nct-iscp) then
!d write (iout,*) 'i.gt.nct-iscp'
iscpend(i,1),*14)
ii=nscp_gr(i)+1
call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,&
- iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),&
+ iatscp_s,iatscp_e,i+iscp,nct_molec(1),nscp_gr(i),iscpstart(i,ii),&
iscpend(i,ii),*14)
endif
enddo ! i
14 continue
#else
iatscp_s=nnt
- iatscp_e=nct-1
- do i=nnt,nct-1
+ iatscp_e=nct_molec(1)-1
+ do i=nnt,nct_molec(1)-1
if (i.lt.nnt+iscp) then
nscp_gr(i)=1
iscpstart(i,1)=i+iscp
- iscpend(i,1)=nct
+ iscpend(i,1)=nct_molec(1)
elseif (i.gt.nct-iscp) then
nscp_gr(i)=1
iscpstart(i,1)=nnt
iscpstart(i,1)=nnt
iscpend(i,1)=i-iscp
iscpstart(i,2)=i+iscp
- iscpend(i,2)=nct
+ iscpend(i,2)=nct_molec(1)
endif
enddo ! i
#endif
endif ! lprint
! Partition local interactions
#ifdef MPI
- call int_bounds(nres-2,loc_start,loc_end)
+ call int_bounds(nres_molec(1)-2,loc_start,loc_end)
loc_start=loc_start+1
loc_end=loc_end+1
- call int_bounds(nres-2,ithet_start,ithet_end)
+ call int_bounds(nres_molec(1)-2,ithet_start,ithet_end)
ithet_start=ithet_start+2
ithet_end=ithet_end+2
- call int_bounds(nct-nnt-2,iturn3_start,iturn3_end)
+ call int_bounds(nct_molec(1)-nnt-2,iturn3_start,iturn3_end)
iturn3_start=iturn3_start+nnt
iphi_start=iturn3_start+2
iturn3_end=iturn3_end+nnt
iphi_end=iturn3_end+2
iturn3_start=iturn3_start-1
iturn3_end=iturn3_end-1
- call int_bounds(nres-3,itau_start,itau_end)
+ call int_bounds(nres_molec(1)-3,itau_start,itau_end)
itau_start=itau_start+3
itau_end=itau_end+3
- call int_bounds(nres-3,iphi1_start,iphi1_end)
+ call int_bounds(nres_molec(1)-3,iphi1_start,iphi1_end)
iphi1_start=iphi1_start+3
iphi1_end=iphi1_end+3
- call int_bounds(nct-nnt-3,iturn4_start,iturn4_end)
+ call int_bounds(nct_molec(1)-nnt-3,iturn4_start,iturn4_end)
iturn4_start=iturn4_start+nnt
iphid_start=iturn4_start+2
iturn4_end=iturn4_end+nnt
iphid_end=iturn4_end+2
iturn4_start=iturn4_start-1
iturn4_end=iturn4_end-1
- call int_bounds(nres-2,ibond_start,ibond_end)
+! print *,"TUTUTU",nres_molec(1),nres
+ call int_bounds(nres_molec(1)-2,ibond_start,ibond_end)
ibond_start=ibond_start+1
ibond_end=ibond_end+1
- call int_bounds(nct-nnt,ibondp_start,ibondp_end)
+ print *,ibond_start,ibond_end
+ call int_bounds(nct_molec(1)-nnt,ibondp_start,ibondp_end)
ibondp_start=ibondp_start+nnt
ibondp_end=ibondp_end+nnt
- call int_bounds1(nres-1,ivec_start,ivec_end)
+ call int_bounds1(nres_molec(1)-1,ivec_start,ivec_end)
! print *,"Processor",myrank,fg_rank,fg_rank1,
! & " ivec_start",ivec_start," ivec_end",ivec_end
iset_start=loc_start+2
iset_end=loc_end+2
- call int_bounds(nres,ilip_start,ilip_end)
+ call int_bounds(nres_molec(1),ilip_start,ilip_end)
ilip_start=ilip_start
ilip_end=ilip_end
- call int_bounds(nres-1,itube_start,itube_end)
+ call int_bounds(nres_molec(1)-1,itube_start,itube_end)
itube_start=itube_start
itube_end=itube_end
if (ndih_constr.eq.0) then
endif
#else
loc_start=2
- loc_end=nres-1
+ loc_end=nres_molec(1)-1
ithet_start=3
- ithet_end=nres
+ ithet_end=nres_molec(1)
iturn3_start=nnt
- iturn3_end=nct-3
+ iturn3_end=nct_molec(1)-3
iturn4_start=nnt
- iturn4_end=nct-4
+ iturn4_end=nct_molec(1)-4
iphi_start=nnt+3
- iphi_end=nct
+ iphi_end=nct_molec(1)
iphi1_start=4
- iphi1_end=nres
+ iphi1_end=nres_molec(1)
idihconstr_start=1
idihconstr_end=ndih_constr
ithetaconstr_start=1
iphid_start=iphi_start
iphid_end=iphi_end-1
itau_start=4
- itau_end=nres
+ itau_end=nres_molec(1)
ibond_start=2
- ibond_end=nres-1
+ ibond_end=nres_molec(1)-1
ibondp_start=nnt
- ibondp_end=nct-1
+ ibondp_end=nct_molec(1)-1
ivec_start=1
- ivec_end=nres-1
+ ivec_end=nres_molec(1)-1
iset_start=3
- iset_end=nres+1
+ iset_end=nres_molec(1)+1
iint_start=2
- iint_end=nres-1
+ iint_end=nres_molec(1)-1
ilip_start=1
- ilip_end=nres
+ ilip_end=nres_molec(1)
itube_start=1
- itube_end=nres
+ itube_end=nres_molec(1)
#endif
!el common /przechowalnia/
! deallocate(iturn3_start_all)
nside=0
do i=2,nres-1
#ifdef WHAM_RUN
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
#else
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
#endif
nside=nside+1
ialph(i,1)=nvar+nside
!-----------------------------------------------------------------------------
! rescode.f
!-----------------------------------------------------------------------------
- integer function rescode(iseq,nam,itype)
+ integer function rescode(iseq,nam,itype,molecule)
use io_base, only: ucase
! implicit real*8 (a-h,o-z)
! include 'COMMON.IOUNITS'
character(len=3) :: nam !,ucase
integer :: iseq,itype,i
-
+ integer :: molecule
+ print *,molecule,nam
+ if (molecule.eq.1) then
if (itype.eq.0) then
- do i=-ntyp1,ntyp1
- if (ucase(nam).eq.restyp(i)) then
+ do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
+ if (ucase(nam).eq.restyp(i,molecule)) then
rescode=i
return
endif
else
- do i=-ntyp1,ntyp1
+ do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
if (nam(1:1).eq.onelet(i)) then
rescode=i
return
enddo
endif
+ else if (molecule.eq.2) then
+ do i=1,ntyp1_molec(molecule)
+ print *,nam(1:1),restyp(i,molecule)(1:1)
+ if (nam(1:1).eq.restyp(i,molecule)(1:1)) then
+ rescode=i
+ return
+ endif
+ enddo
+ else if (molecule.eq.3) then
+ write(iout,*) "SUGAR not yet implemented"
+ stop
+ else if (molecule.eq.4) then
+ write(iout,*) "Explicit LIPID not yet implemented"
+ stop
+ else if (molecule.eq.5) then
+ do i=1,ntyp1_molec(molecule)
+ print *,i,restyp(i,molecule)
+ if (ucase(nam).eq.restyp(i,molecule)) then
+ rescode=i
+ return
+ endif
+ enddo
+ else
+ write(iout,*) "molecule not defined"
+ endif
write (iout,10) iseq,nam
stop
10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
end function rescode
+ integer function sugarcode(sugar,ires)
+ character sugar
+ integer ires
+ if (sugar.eq.'D') then
+ sugarcode=1
+ else if (sugar.eq.' ') then
+ sugarcode=2
+ else
+ write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
+ stop
+ endif
+ return
+ end function sugarcode
+
!-----------------------------------------------------------------------------
! timing.F
!-----------------------------------------------------------------------------
logical :: minim,refstr,pdbref,overlapsc,&
energy_dec,sideadd,lsecondary,read_cart,unres_pdb,&
vdisulf,searchsc,lmuca,dccart,extconf,out1file,&
- gnorm_check,gradout,split_ene,with_theta_constr
+ gnorm_check,gradout,split_ene,with_theta_constr,protein,ions,nucleic
#ifdef CLUSTER
integer :: iopt,nend,nstart,outpdb,outmol2 !cluster
logical :: punch_dist,print_dist,lside,lprint_cart,lprint_int,&
real(kind=8),dimension(2,2) :: app,bpp,ael6,ael3
integer :: expon,expon2, nnt,nct,itypro
integer,dimension(:,:),allocatable :: istart,iend !(maxres,maxint_gr)
- integer,dimension(:),allocatable :: nint_gr,itype,itel,&
+ integer,dimension(:),allocatable :: nint_gr,itel,&
ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr !(maxres)
+ integer,dimension(:),allocatable :: istype
+ integer,dimension(:,:),allocatable :: itype ! now itype has more molecule types
integer,dimension(:,:),allocatable :: iscpstart,iscpend !(maxres,maxint_gr)
integer :: iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,&
iatel_e_vdw,iatscp_s,iatscp_e,ispp,iscp
real(kind=8),dimension(:,:),allocatable :: xloc,xrot !(3,maxres)
real(kind=8),dimension(:),allocatable :: dc_work !(MAXRES6)
integer :: nres,nres0
+ integer, dimension(5) :: nres_molec
! common /rotmat/
real(kind=8),dimension(:,:,:),allocatable :: prod,rt !(3,3,maxres)
! common /refstruct/
real(kind=8) :: buflipbot, bufliptop,bordlipbot,bordliptop, &
lipbufthick,lipthick
integer,dimension(:,:),allocatable :: tabperm !(maxperm,maxsym)
+ integer molec
! common /from_zscore/ in module.compare
!-----------------------------------------------------------------------------
! common.geo
!-----------------------------------------------------------------------------
! Number of AA types (at present only natural AA's will be handled
integer,parameter :: ntyp=24,ntyp1=ntyp+1
+ integer,dimension(5) :: ntyp_molec=(/24,5,0,0,5/),ntyp1_molec=(/25,6,0,0,6/)
+ integer,parameter ::maxmolec=5
+
!-----------------------------------------------------------------------------
! common.names
! common /names/
!-----------------------------------------------------------------------------
! block data nazwy
!el allocate(restyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
- character(len=3),dimension(-ntyp1:ntyp1) :: restyp = &
+ character(len=3),dimension(-ntyp1:ntyp1,maxmolec) :: restyp = &
(/'DD ','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS',&
'DGL','DSG','DGN','DSN','DTH',&
'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',&
'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',&
'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',&
- 'AIB','ABU','D '/)
+ 'AIB','ABU','D ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',&
+ 'A ','G ','C ','T ','U ','X ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',&
+ 'NA+','MG2','K+ ','CA2','CL-',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',&
+ ' ',' ',' ',' ',' ',' ',' ',' ',' ',' '&
+ /)
!el allocate(onelet(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
character(len=1),dimension(-ntyp1:ntyp1) :: onelet = &
(/'z','z','z','z','z','p','k','r','h','d','e','n','q','s',&
't','g','a','y','w','v','l','i','f','m','c','x',&
'C','M','F','I','L','V','W','Y','A','G','T',&
'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/)
+! character(len=1),dimension(ntyp1_nucl) :: restyp_nucl = &
+! (/'A','G','C','T','U','X'/)
+
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! Number of energy components
integer :: nprint_ene = 21
integer,dimension(n_ene) :: print_order = &
(/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25/)
+
+
!#endif
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! Calculate the bond-stretching energy
!
call ebond(estr)
+ print *,"EBOND",estr
! write(iout,*) "in etotal afer ebond",ipot
!
! allocate(gacont(3,nres/4,iatsc_s:iatsc_e)) !(3,maxconts,maxres)
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
evdw=evdw+evdwij
! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
! endif
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
!d & om1,om2,om12,1.0D0/dsqrt(rrij),
!el ind=0
do i=iatsc_s,iatsc_e
!C print *,"I am in EVDW",i
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
! if (i.ne.47) cycle
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
enddo! k
ELSE
!el ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! if (j.ne.78) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,&
! 1.0d0/vbld(j+nres) !d
-! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
if (rij_shift.le.0.0D0) then
evdw=1.0D20
!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
return
endif
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa!(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
- restyp(itypi),i,restyp(itypj),j, &
+ restyp(itypi,1),i,restyp(itypj,1),j, &
epsi,sigm,chi1,chi2,chip1,chip2, &
eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
! if (icall.eq.0) lprn=.true.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
- restyp(itypi),i,restyp(itypj),j,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,&
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
eello_turn4=0.0d0
!el ind=0
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_conti=0
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
!el ind=ind+1
iteli=itel(i)
itelj=itel(j)
endif
! if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
if (i.gt. nnt+2 .and. i.lt.nct+2) then
- iti = itortyp(itype(i-2))
+ iti = itortyp(itype(i-2,1))
else
iti=ntortyp+1
endif
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
- iti1 = itortyp(itype(i-1))
+ iti1 = itortyp(itype(i-1,1))
else
iti1=ntortyp+1
endif
-! print *,iti,i,"iti",iti1,itype(i-1),itype(i-2)
+! print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1)
!d write (iout,*) '*******i',i,' iti1',iti
!d write (iout,*) 'b1',b1(:,iti)
!d write (iout,*) 'b2',b2(:,iti)
enddo
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
- if (itype(i-1).le.ntyp) then
- iti1 = itortyp(itype(i-1))
+ if (itype(i-1,1).le.ntyp) then
+ iti1 = itortyp(itype(i-1,1))
else
iti1=ntortyp+1
endif
#endif
#endif
!d do i=1,nres
-!d iti = itortyp(itype(i))
+!d iti = itortyp(itype(i,1))
!d write (iout,*) i
!d do j=1,2
!d write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)')
! print *,"before iturn3 loop"
do i=iturn3_start,iturn3_end
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
- .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
- .or. itype(i+3).eq.ntyp1 &
- .or. itype(i+4).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+3,1).eq.ntyp1 &
+ .or. itype(i+4,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_conti=num_cont_hb(i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
- if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
+ if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
enddo ! i
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
-! write (iout,*) i,j,itype(i),itype(j)
- if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+! write (iout,*) i,j,itype(i,1),itype(j,1)
+ if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
a32=a32*fac
a33=a33*fac
!d write (iout,'(4i5,4f10.5)')
-!d & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
!d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
!d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
!d & uy(:,j),uz(:,j)
a_temp(1,2)=a23
a_temp(2,1)=a32
a_temp(2,2)=a33
- iti1=itortyp(itype(i+1))
- iti2=itortyp(itype(i+2))
- iti3=itortyp(itype(i+3))
+ iti1=itortyp(itype(i+1,1))
+ iti2=itortyp(itype(i+2,1))
+ iti3=itortyp(itype(i+3,1))
! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
call transpose2(EUg(1,1,i+1),e1t(1,1))
call transpose2(Eug(1,1,i+2),e2t(1,1))
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- if (itype(j).eq.ntyp1) cycle
- itypj=iabs(itype(j))
+ if (itype(j,1).eq.ntyp1) cycle
+ itypj=iabs(itype(j,1))
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
! yj=c(2,nres+j)-yi
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
! 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
if (.not.dyn_ss .and. i.le.nss) then
! 15/02/13 CC dynamic SSbond - additional check
- if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. &
- iabs(itype(jjj)).eq.1) then
+ if (ii.gt.nres .and. iabs(itype(iii,1)).eq.1 .and. &
+ iabs(itype(jjj,1)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
!d write (iout,*) "eij",eij
deltat1,deltat2,deltat12,ed,pom1,pom2,eom1,eom2,eom12,&
cosphi,ggk
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
dzi=dc_norm(3,nres+i)
! dsci_inv=dsc_inv(itypi)
dsci_inv=vbld_inv(nres+i)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(nres+j)
xj=c(1,nres+j)-xi
! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
do i=ibondp_start,ibondp_end
- if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
- if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+ if (itype(i-1,1).eq.ntyp1 .and. itype(i,1).eq.ntyp1) cycle
+ if (itype(i-1,1).eq.ntyp1 .or. itype(i,1).eq.ntyp1) then
!C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
!C do j=1,3
!C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
! endif
enddo
estr=0.5d0*AKP*estr+estr1
+ print *,"estr_bb",estr,AKP
!
! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
!
do i=ibond_start,ibond_end
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
+ if (iti.eq.0) print *,"WARNING WRONG SETTTING",i
if (iti.ne.10 .and. iti.ne.ntyp1) then
nbi=nbondterm(iti)
if (nbi.eq.1) then
"estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
+ print *,"estr_sc",estr
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
enddo
usumsqder=usumsqder+ud(j)*uprod2
enddo
estr=estr+uprod/usum
+ print *,"estr_sc",estr,i
+
if (energy_dec) write (iout,*) &
"estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
AKSC(1,iti),AKSC(1,iti)*diff*diff
etheta=0.0D0
! write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+ if (itype(i-1,1).eq.ntyp1) cycle
! Zero the energy function and its derivative at 0 or pi.
call splinthet(theta(i),0.5d0*delta,ss,ssd)
- it=itype(i-1)
- ichir1=isign(1,itype(i-2))
- ichir2=isign(1,itype(i))
- if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
- if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
- if (itype(i-1).eq.10) then
- itype1=isign(10,itype(i-2))
- ichir11=isign(1,itype(i-2))
- ichir12=isign(1,itype(i-2))
- itype2=isign(10,itype(i))
- ichir21=isign(1,itype(i))
- ichir22=isign(1,itype(i))
+ it=itype(i-1,1)
+ ichir1=isign(1,itype(i-2,1))
+ ichir2=isign(1,itype(i,1))
+ if (itype(i-2,1).eq.10) ichir1=isign(1,itype(i-1,1))
+ if (itype(i,1).eq.10) ichir2=isign(1,itype(i-1,1))
+ if (itype(i-1,1).eq.10) then
+ itype1=isign(10,itype(i-2,1))
+ ichir11=isign(1,itype(i-2,1))
+ ichir12=isign(1,itype(i-2,1))
+ itype2=isign(10,itype(i,1))
+ ichir21=isign(1,itype(i,1))
+ ichir22=isign(1,itype(i,1))
endif
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.gt.3 .and. itype(i-2,1).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
y(1)=0.0D0
y(2)=0.0D0
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i,1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
etheta=0.0D0
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
- if (itype(i-2).eq.ntyp1.or.itype(i).eq.ntyp1) cycle
- if (iabs(itype(i+1)).eq.20) iblock=2
- if (iabs(itype(i+1)).ne.20) iblock=1
+ if (itype(i-1,1).eq.ntyp1) cycle
+ if (itype(i-2,1).eq.ntyp1.or.itype(i,1).eq.ntyp1) cycle
+ if (iabs(itype(i+1,1)).eq.20) iblock=2
+ if (iabs(itype(i+1,1)).ne.20) iblock=1
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
theti2=0.5d0*theta(i)
- ityp2=ithetyp((itype(i-1)))
+ ityp2=ithetyp((itype(i-1,1)))
do k=1,nntheterm
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1),1).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
#else
phii=phi(i)
#endif
- ityp1=ithetyp((itype(i-2)))
+ ityp1=ithetyp((itype(i-2,1)))
! propagation of chirality for glycine type
do k=1,nsingle
cosph1(k)=dcos(k*phii)
enddo
else
phii=0.0d0
- ityp1=ithetyp(itype(i-2))
+ ityp1=ithetyp(itype(i-2,1))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i+1,1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
#else
phii1=phi(i+1)
#endif
- ityp3=ithetyp((itype(i)))
+ ityp3=ithetyp((itype(i,1)))
do k=1,nsingle
cosph2(k)=dcos(k*phii1)
sinph2(k)=dsin(k*phii1)
enddo
else
phii1=0.0d0
- ityp3=ithetyp(itype(i))
+ ityp3=ithetyp(itype(i,1))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
escloc=0.0D0
! write (iout,'(a)') 'ESC'
do i=loc_start,loc_end
- it=itype(i)
+ it=itype(i,1)
if (it.eq.ntyp1) cycle
if (it.eq.10) goto 1
nlobit=nlob(iabs(it))
delta=0.02d0*pi
escloc=0.0D0
do i=loc_start,loc_end
- if (itype(i).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1) cycle
costtab(i+1) =dcos(theta(i+1))
sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=iabs(itype(i))
+ it=iabs(itype(i,1))
if (it.eq.10) goto 1
!
! Compute the axes of tghe local cartesian coordinates system; store in
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
do j = 1,3
- z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
+ z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i,1)))
enddo
! write (2,*) "i",i
! write (2,*) "x_prime",(x_prime(j),j=1,3)
! Compute the energy of the ith side cbain
!
! write (2,*) "xx",xx," yy",yy," zz",zz
- it=iabs(itype(i))
+ it=iabs(itype(i,1))
do j = 1,65
x(j) = sc_parmin(j,it)
enddo
!c diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0,dfloat(itype(i,1)))*dsin(alph(2))*dsin(omeg(2))
write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
xx1,yy1,zz1
! & dscp1,dscp2,sumene
! sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
-! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1)
! & ,zz,xx,yy
!#define DEBUG
#ifdef DEBUG
!
! Compute the gradient of esc
!
-! zz=zz*dsign(1.0,dfloat(itype(i)))
+! zz=zz*dsign(1.0,dfloat(itype(i,1)))
pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
+(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) &
+(pom1+pom2)*pom_dx
#ifdef DEBUG
- write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i)
+ write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i,1)
#endif
!
sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
+(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) &
+(pom1-pom2)*pom_dy
#ifdef DEBUG
- write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i)
+ write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i,1)
#endif
!
de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy &
+x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) &
+ ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6)
#ifdef DEBUG
- write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i)
+ write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i,1)
#endif
!
de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) &
-0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) &
+pom1*pom_dt1+pom2*pom_dt2
#ifdef DEBUG
- write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
+ write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i,1)
#endif
!
!
dZZ_Ci(k)=0.0d0
do j=1,3
dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) &
- *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres)
dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) &
- *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres)
enddo
dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
etors=0.0D0
do i=iphi_start,iphi_end
etors_ii=0.0D0
- if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 &
- .or. itype(i).eq.ntyp1) cycle
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
+ if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1) cycle
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
phii=phi(i)
gloci=0.0D0
! Proline-Proline pair is a special case...
'etor',i,etors_ii
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
- restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,&
+ restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,&
(v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
! write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
! lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
- .or. itype(i-3).eq.ntyp1 &
- .or. itype(i).eq.ntyp1) cycle
+ if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i-3,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1) cycle
etors_ii=0.0D0
- if (iabs(itype(i)).eq.20) then
+ if (iabs(itype(i,1)).eq.20) then
iblock=2
else
iblock=1
endif
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
phii=phi(i)
gloci=0.0D0
! Regular cosine and sine terms
'etor',i,etors_ii-v0(itori,itori1,iblock)
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
- restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,&
+ restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,&
(v1(j,itori,itori1,iblock),j=1,6),&
(v2(j,itori,itori1,iblock),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
! write(iout,*) "a tu??"
do i=iphid_start,iphid_end
etors_d_ii=0.0D0
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
- .or. itype(i-3).eq.ntyp1 &
- .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
- itori2=itortyp(itype(i))
+ if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i-3,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
+ itori2=itortyp(itype(i,1))
phii=phi(i)
phii1=phi(i+1)
gloci1=0.0D0
gloci2=0.0D0
iblock=1
- if (iabs(itype(i+1)).eq.20) iblock=2
+ if (iabs(itype(i+1,1)).eq.20) iblock=2
! Regular cosine and sine terms
do j=1,ntermd_1(itori,itori1,itori2,iblock)
! write (iout,*) "EBACK_SC_COR",itau_start,itau_end
esccor=0.0D0
do i=itau_start,itau_end
- if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
+ if ((itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1)) cycle
esccor_ii=0.0D0
- isccori=isccortyp(itype(i-2))
- isccori1=isccortyp(itype(i-1))
+ isccori=isccortyp(itype(i-2,1))
+ isccori1=isccortyp(itype(i-1,1))
! write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
phii=phi(i)
! 2 = Ca...Ca...Ca...SC
! 3 = SC...Ca...Ca...SCi
gloci=0.0D0
- if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. &
- (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. &
- (itype(i-1).eq.ntyp1))) &
- .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) &
- .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) &
- .or.(itype(i).eq.ntyp1))) &
- .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. &
- (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. &
- (itype(i-3).eq.ntyp1)))) cycle
- if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
- if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) &
+ if (((intertyp.eq.3).and.((itype(i-2,1).eq.10).or. &
+ (itype(i-1,1).eq.10).or.(itype(i-2,1).eq.ntyp1).or. &
+ (itype(i-1,1).eq.ntyp1))) &
+ .or. ((intertyp.eq.1).and.((itype(i-2,1).eq.10) &
+ .or.(itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1) &
+ .or.(itype(i,1).eq.ntyp1))) &
+ .or.((intertyp.eq.2).and.((itype(i-1,1).eq.10).or. &
+ (itype(i-1,1).eq.ntyp1).or.(itype(i-2,1).eq.ntyp1).or. &
+ (itype(i-3,1).eq.ntyp1)))) cycle
+ if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1,1).eq.ntyp1)) cycle
+ if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres,1).eq.ntyp1)) &
cycle
do j=1,nterm_sccor(isccori,isccori1)
v1ij=v1sccor(j,intertyp,isccori,isccori1)
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
- restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,&
+ restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,isccori,isccori1,&
(v1sccor(j,intertyp,isccori,isccori1),j=1,6),&
(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
allocate(dipderx(3,5,4,maxconts,nres))
!
- iti1 = itortyp(itype(i+1))
+ iti1 = itortyp(itype(i+1,1))
if (j.lt.nres-1) then
- itj1 = itortyp(itype(j+1))
+ itj1 = itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
if (l.eq.j+1) then
! parallel orientation of the two CA-CA-CA frames.
if (i.gt.1) then
- iti=itortyp(itype(i))
+ iti=itortyp(itype(i,1))
else
iti=ntortyp+1
endif
- itk1=itortyp(itype(k+1))
- itj=itortyp(itype(j))
+ itk1=itortyp(itype(k+1,1))
+ itj=itortyp(itype(j,1))
if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
+ itl1=itortyp(itype(l+1,1))
else
itl1=ntortyp+1
endif
else
! Antiparallel orientation of the two CA-CA-CA frames.
if (i.gt.1) then
- iti=itortyp(itype(i))
+ iti=itortyp(itype(i,1))
else
iti=ntortyp+1
endif
- itk1=itortyp(itype(k+1))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+ itk1=itortyp(itype(k+1,1))
+ itl=itortyp(itype(l,1))
+ itj=itortyp(itype(j,1))
if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
+ itj1=itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
!d write (iout,*)
!d & 'EELLO5: Contacts have occurred for peptide groups',i,j,
!d & ' and',k,l
- itk=itortyp(itype(k))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+ itk=itortyp(itype(k,1))
+ itl=itortyp(itype(l,1))
+ itj=itortyp(itype(j,1))
eello5_1=0.0d0
eello5_2=0.0d0
eello5_3=0.0d0
! i i C
! C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- itk=itortyp(itype(k))
+ itk=itortyp(itype(k,1))
s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
!
! 4/7/01 AL Component s1 was removed, because it pertains to the respective
! energy moment and not to the cluster cumulant.
- iti=itortyp(itype(i))
+ iti=itortyp(itype(i,1))
if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
+ itj1=itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
- itk=itortyp(itype(k))
- itk1=itortyp(itype(k+1))
+ itk=itortyp(itype(k,1))
+ itk1=itortyp(itype(k+1,1))
if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
+ itl1=itortyp(itype(l+1,1))
else
itl1=ntortyp+1
endif
! 4/7/01 AL Component s1 was removed, because it pertains to the respective
! energy moment and not to the cluster cumulant.
!d write (2,*) 'eello_graph4: wturn6',wturn6
- iti=itortyp(itype(i))
- itj=itortyp(itype(j))
+ iti=itortyp(itype(i,1))
+ itj=itortyp(itype(j,1))
if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
+ itj1=itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
- itk=itortyp(itype(k))
+ itk=itortyp(itype(k,1))
if (k.lt.nres-1) then
- itk1=itortyp(itype(k+1))
+ itk1=itortyp(itype(k+1,1))
else
itk1=ntortyp+1
endif
- itl=itortyp(itype(l))
+ itl=itortyp(itype(l,1))
if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
+ itl1=itortyp(itype(l+1,1))
else
itl1=ntortyp+1
endif
j=i+4
k=i+1
l=i+3
- iti=itortyp(itype(i))
- itk=itortyp(itype(k))
- itk1=itortyp(itype(k+1))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+ iti=itortyp(itype(i,1))
+ itk=itortyp(itype(k,1))
+ itk1=itortyp(itype(k+1,1))
+ itl=itortyp(itype(l,1))
+ itj=itortyp(itype(j,1))
!d write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj
!d write (2,*) 'i',i,' k',k,' j',j,' l',l
!d if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
! Derivatives in alpha and omega:
!
do i=2,nres-1
-! dsci=dsc(itype(i))
+! dsci=dsc(itype(i,1))
dsci=vbld(i+nres)
#ifdef OSF
alphi=alph(i)
! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
! endif
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
!d & om1,om2,om12,1.0D0/dsqrt(rrij),
! endif
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
!d & om1,om2,om12,1.0D0/dsqrt(rrij),
! if (icall.eq.0) lprn=.false.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
ELSE
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
! & 1.0d0/vbld(j+nres)
-! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
if (rij_shift.le.0.0D0) then
evdw=1.0D20
!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
return
endif
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
- restyp(itypi),i,restyp(itypj),j,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
! if (icall.eq.0) lprn=.false.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
! 'evdw',i,j,evdwij,' ss'
ELSE
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
! & 1.0d0/vbld(j+nres)
-! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
if (rij_shift.le.0.0D0) then
evdw=1.0D20
!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
return
endif
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
- restyp(itypi),i,restyp(itypj),j,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
! if (icall.eq.0) lprn=.true.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
- restyp(itypi),i,restyp(itypj),j,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,&
! if (icall.eq.0) lprn=.true.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
- restyp(itypi),i,restyp(itypj),j,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,&
! Loop over i,i+2 and i,i+3 pairs of the peptide groups
!
do i=iturn3_start,iturn3_end
- if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 &
- .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
- .or. itype(i+3).eq.ntyp1 &
- .or. itype(i+4).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+3,1).eq.ntyp1 &
+ .or. itype(i+4,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=num_cont_hb(i)
call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
- if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
+ if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
enddo ! i
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
call eelecij_scale(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
a32=a32*fac
a33=a33*fac
!d write (iout,'(4i5,4f10.5)')
-!d & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
!d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
!d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
!d & uy(:,j),uz(:,j)
! & " iatel_e_vdw",iatel_e_vdw
call flush(iout)
do i=iatel_s_vdw,iatel_e_vdw
- if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
! & ' ielend',ielend_vdw(i)
call flush(iout)
do j=ielstart_vdw(i),ielend_vdw(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
!el ind=ind+1
iteli=itel(i)
itelj=itel(j)
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
enddo
if (n.le.nphi+ntheta) goto 10
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
galphai=0.0D0
gomegai=0.0D0
do k=1,3
do j=1,3
dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/&
vbld(i-1)
- if (itype(i-1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
+ if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/&
vbld(i)
- if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
+ if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
enddo
enddo
#if defined(MPI) && defined(PARINTDER)
#else
do i=3,nres
#endif
- if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
+ if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1)) then
cost1=dcos(omicron(1,i))
sint1=sqrt(1-cost1*cost1)
cost2=dcos(omicron(2,i))
dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) &
+cost2*(-dc_norm(j,i-1+nres)))/ &
vbld(i-1+nres)
-! write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
+! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres)
domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
enddo
endif
#else
do i=4,nres
#endif
-! if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle
+! if (itype(i-1,1).eq.21 .or. itype(i-2,1).eq.21 ) cycle
! the conventional case
sint=dsin(theta(i))
sint1=dsin(theta(i-1))
ctgt=cost/sint
ctgt1=cost1/sint1
cosg_inv=1.0d0/cosg
- if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+ if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
-(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
! Obtaining the gamma derivatives from cosine derivative
else
do j=1,3
- if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+ if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
dc_norm(j,i-3))/vbld(i-2)
do i=3,nres
!elwrite(iout,*) " vecpr",i,nres
#endif
- if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
-! if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
-! & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
+ if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
+! if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10).or.
+! & (itype(i-1,1).eq.ntyp1).or.(itype(i,1).eq.ntyp1)) cycle
!c dtauangle(j,intertyp,dervityp,residue number)
!c INTERTYP=1 SC...Ca...Ca..Ca
! the conventional case
#else
do i=4,nres
#endif
- if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. &
- (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
+ if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+ (itype(i-2,1).eq.ntyp1).or.(itype(i-3,1).eq.ntyp1)) cycle
! the conventional case
sint=dsin(omicron(1,i))
sint1=dsin(theta(i-1))
do i=3,nres
#endif
! the conventional case
- if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. &
- (itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
+ if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+ (itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
sint=dsin(omicron(1,i))
sint1=dsin(omicron(2,i-1))
sing=dsin(tauangle(3,i))
#else
do i=2,nres-1
#endif
- if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
fac6=fac5/vbld(i)
fac7=fac5*fac5
write (iout,*) &
"Analytical (upper) and numerical (lower) gradient of alpha"
do i=2,nres-1
- if(itype(i).ne.10) then
+ if(itype(i,1).ne.10) then
do j=1,3
dcji=dc(j,i-1)
dc(j,i-1)=dcji+aincr
write (iout,*) &
"Analytical (upper) and numerical (lower) gradient of omega"
do i=2,nres-1
- if(itype(i).ne.10) then
+ if(itype(i,1).ne.10) then
do j=1,3
dcji=dc(j,i-1)
dc(j,i-1)=dcji+aincr
(cref(3,jl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,jl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
(cref(3,jl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,jl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
dqwol(k,jl)=dqwol(k,jl)-ddqij
enddo
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
dqwol(k,il)=dqwol(k,il)+ddqij
dqwol(k,jl)=dqwol(k,jl)-ddqij
enddo
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
!el allocate(dyn_ssbond_ij(iatsc_s:iatsc_e,nres))
!el allocate(dyn_ssbond_ij(0:nres+4,nres))
- itypi=itype(i)
+ itypi=itype(i,1)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
- itypj=itype(j)
+ itypj=itype(j,1)
xj=c(1,nres+j)-c(1,nres+i)
yj=c(2,nres+j)-c(2,nres+i)
zj=c(3,nres+j)-c(3,nres+i)
j=resj
k=resk
!C write(iout,*) resi,resj,resk
- itypi=itype(i)
+ itypi=itype(i,1)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
- itypj=itype(j)
+ itypj=itype(j,1)
xj=c(1,nres+j)
yj=c(2,nres+j)
zj=c(3,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
dscj_inv=vbld_inv(j+nres)
- itypk=itype(k)
+ itypk=itype(k,1)
xk=c(1,nres+k)
yk=c(2,nres+k)
zk=c(3,nres+k)
! print *, "I am in eliptran"
do i=ilip_start,ilip_end
!C do i=1,1
- if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))&
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1).or.(i.eq.nres))&
cycle
positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
enddo
! here starts the side chain transfer
do i=ilip_start,ilip_end
- if (itype(i).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1) cycle
positi=(mod(c(3,i+nres),boxzsize))
if (positi.le.0) positi=positi+boxzsize
!C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
!C lipbufthick is thickenes of lipid buffore
sslip=sscalelip(fracinbuf)
ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
- eliptran=eliptran+sslip*liptranene(itype(i))
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
gliptranx(3,i)=gliptranx(3,i) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
gliptranc(3,i-1)= gliptranc(3,i-1) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
!C print *,"doing sccale for lower part"
elseif (positi.gt.bufliptop) then
fracinbuf=1.0d0- &
((bordliptop-positi)/lipbufthick)
sslip=sscalelip(fracinbuf)
ssgradlip=sscagradlip(fracinbuf)/lipbufthick
- eliptran=eliptran+sslip*liptranene(itype(i))
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
gliptranx(3,i)=gliptranx(3,i) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
gliptranc(3,i-1)= gliptranc(3,i-1) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
!C print *, "doing sscalefor top part",sslip,fracinbuf
else
- eliptran=eliptran+liptranene(itype(i))
+ eliptran=eliptran+liptranene(itype(i,1))
!C print *,"I am in true lipid"
endif
endif ! if in lipid or buffor
!C for UNRES
do i=itube_start,itube_end
!C lets ommit dummy atoms for now
- if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
!C now calculate distance from center of tube and direction vectors
xmin=boxxsize
ymin=boxysize
do i=itube_start,itube_end
!C Lets not jump over memory as we use many times iti
- iti=itype(i)
+ iti=itype(i,1)
!C lets ommit dummy atoms for now
if ((iti.eq.ntyp1) &
!C in UNRES uncomment the line below as GLY has no side-chain...
do i=itube_start,itube_end
!C lets ommit dummy atoms for now
- if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
!C now calculate distance from center of tube and direction vectors
!C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
!C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
!C lipbufthick is thickenes of lipid buffore
sstube=sscalelip(fracinbuf)
ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-!C print *,ssgradtube, sstube,tubetranene(itype(i))
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
enetube(i)=enetube(i)+sstube*tubetranenepep
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *,"doing sccale for lower part"
elseif (positi.gt.buftubetop) then
fracinbuf=1.0d0- &
ssgradtube=sscagradlip(fracinbuf)/tubebufthick
enetube(i)=enetube(i)+sstube*tubetranenepep
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *, "doing sscalefor top part",sslip,fracinbuf
else
sstube=1.0d0
!C print *,gg_tube(1,0),"TU"
do i=itube_start,itube_end
!C Lets not jump over memory as we use many times iti
- iti=itype(i)
+ iti=itype(i,1)
!C lets ommit dummy atoms for now
if ((iti.eq.ntyp1) &
!!C in UNRES uncomment the line below as GLY has no side-chain...
!C lipbufthick is thickenes of lipid buffore
sstube=sscalelip(fracinbuf)
ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-!C print *,ssgradtube, sstube,tubetranene(itype(i))
- enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *,"doing sccale for lower part"
elseif (positi.gt.buftubetop) then
fracinbuf=1.0d0- &
sstube=sscalelip(fracinbuf)
ssgradtube=sscagradlip(fracinbuf)/tubebufthick
- enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *, "doing sscalefor top part",sslip,fracinbuf
else
sstube=1.0d0
ssgradtube=0.0d0
- enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
!C print *,"I am in true lipid"
endif
else
!C for UNRES
do i=itube_start,itube_end
!C lets ommit dummy atoms for now
- if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
!C now calculate distance from center of tube and direction vectors
xmin=boxxsize
ymin=boxysize
do i=itube_start,itube_end
enecavtube(i)=0.0d0
!C Lets not jump over memory as we use many times iti
- iti=itype(i)
+ iti=itype(i,1)
!C lets ommit dummy atoms for now
if ((iti.eq.ntyp1) &
!C in UNRES uncomment the line below as GLY has no side-chain...
enddo
do i=ivec_start,ivec_end
!C do i=1,nres-1
-!C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+!C if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
ishield_list(i)=0
- if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
!Cif there two consequtive dummy atoms there is no peptide group between them
!C the line below has to be changed for FGPROC>1
VolumeTotal=0.0
do k=1,nres
- if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ if ((itype(k,1).eq.ntyp1).or.(itype(k,1).eq.10)) cycle
dist_pep_side=0.0
dist_side_calf=0.0
do j=1,3
enddo
endif
!C this is what is now we have the distance scaling now volume...
- short=short_r_sidechain(itype(k))
- long=long_r_sidechain(itype(k))
+ short=short_r_sidechain(itype(k,1))
+ long=long_r_sidechain(itype(k,1))
costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
sinthet=short/dist_pep_side*costhet
!C now costhet_grad
enddo
fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
-!C write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
+!C write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i)
enddo
return
end subroutine set_shield_fac2
- module geometry
+ module geometry
!-----------------------------------------------------------------------------
use io_units
use names
be1=rad2deg*beta(nres+i,i,nres2+2,i+1)
alfai=0.0D0
if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
- write (iout,1212) restyp(itype(i)),i,dist(i-1,i),&
+ write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),&
alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
enddo
1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
real(kind=8) :: alphi,omegi,theta2
real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
real(kind=8) :: xp,yp,zp,cost2,sint2,rj
-! dsci=dsc(itype(i))
-! dsci_inv=dsc_inv(itype(i))
+! dsci=dsc(itype(i,1))
+! dsci_inv=dsc_inv(itype(i,1))
dsci=vbld(i+nres)
dsci_inv=vbld_inv(i+nres)
#ifdef OSF
xx(1)= xp*cost2+yp*sint2
xx(2)=-xp*sint2+yp*cost2
xx(3)= zp
-!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
+!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i,1)),i,
!d & xp,yp,zp,(xx(k),k=1,3)
do j=1,3
xloc(j,i)=xx(j)
be=0.0D0
if (i.gt.2) then
if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
- if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
+ if ((itype(i,1).ne.10).and.(itype(i-1,1).ne.10)) then
tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
endif
- if (itype(i-1).ne.10) then
+ if (itype(i-1,1).ne.10) then
tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
omicron(1,i)=alpha(i-2,i-1,i-1+nres)
omicron(2,i)=alpha(i-1+nres,i-1,i)
endif
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
endif
endif
vbld(i)=dist(i-1,i)
vbld_inv(i)=1.0d0/vbld(i)
vbld(nres+i)=dist(nres+i,i)
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
vbld_inv(nres+i)=1.0d0/vbld(nres+i)
else
vbld_inv(nres+i)=0.0d0
enddo
if (lprn) then
do i=2,nres
- write (iout,1212) restyp(itype(i)),i,vbld(i),&
+ write (iout,1212) restyp(itype(i,1),1),i,vbld(i),&
rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),&
rad2deg*alph(i),rad2deg*omeg(i)
enddo
print *,'dv=',dv
do 10 it=1,1
if (it.eq.10) goto 10
- open (20,file=restyp(it)//'_distr.sdc',status='unknown')
+ open (20,file=restyp(it,1)//'_distr.sdc',status='unknown')
call gen_side(it,90.0D0 * deg2rad,al,om,fail)
close (20)
goto 10
- open (20,file=restyp(it)//'_distr1.sdc',status='unknown')
+ open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown')
do i=0,90
do j=0,72
prob(j,i)=0.0D0
maxsi=100
!d write (iout,*) 'Gen_Rand_conf: nstart=',nstart
if (nstart.lt.5) then
- it1=iabs(itype(2))
- phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
+ it1=iabs(itype(2,1))
+ phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1)))
! write(iout,*)'phi(4)=',rad2deg*phi(4)
- if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
+ if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4))
! write(iout,*)'theta(3)=',rad2deg*theta(3)
if (it1.ne.10) then
nsi=0
endif
return 1
endif
- it1=iabs(itype(i-1))
- it2=iabs(itype(i-2))
- it=iabs(itype(i))
+ it1=iabs(itype(i-1,1))
+ it2=iabs(itype(i-2,1))
+ it=iabs(itype(i,1))
! print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
! & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
phi(i+1)=gen_phi(i+1,it1,it)
nres2=2*nres
data redfac /0.5D0/
overlap=.false.
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
if (iti.gt.ntyp) return
! Check for SC-SC overlaps.
!d print *,'nnt=',nnt,' nct=',nct
do j=nnt,i-1
- itj=iabs(itype(j))
+ itj=iabs(itype(j,1))
if (j.lt.i-1 .or. ipot.ne.4) then
rcomp=sigmaii(iti,itj)
else
c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
enddo
do j=nnt,i-2
- itj=iabs(itype(j))
+ itj=iabs(itype(j,1))
!d print *,'overlap, p-Sc: i=',i,' j=',j,
!d & ' dist=',dist(nres+j,maxres2+1)
if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
do ires=1,ioverlap_last
i=ioverlap(ires)
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
if (iti.ne.10) then
nsi=0
fail=.true.
! print *,'>>overlap_sc nnt=',nnt,' nct=',nct
ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
- itypi1=iabs(itype(i+1))
+ itypi=iabs(itype(i,1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
dscj_inv=dsc_inv(itypj)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
endif
do i=1,nres-1
!in wham do i=1,nres
- iti=itype(i)
+ iti=itype(i,1)
if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
- (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1)) then
+ (iti.ne.ntyp1 .and. itype(i+1,1).ne.ntyp1)) then
write (iout,'(a,i4)') 'Bad Cartesians for residue',i
!test stop
endif
enddo
!el -----
!#ifdef WHAM_RUN
-! if (itype(1).eq.ntyp1) then
+! if (itype(1,1).eq.ntyp1) then
! do j=1,3
! c(j,1)=c(j,2)+(c(j,3)-c(j,4))
! enddo
! endif
-! if (itype(nres).eq.ntyp1) then
+! if (itype(nres,1).eq.ntyp1) then
! do j=1,3
! c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
! enddo
! endif
!#endif
! if (unres_pdb) then
-! if (itype(1).eq.21) then
+! if (itype(1,1).eq.21) then
! theta(3)=90.0d0*deg2rad
! phi(4)=180.0d0*deg2rad
! vbld(2)=3.8d0
! vbld_inv(2)=1.0d0/vbld(2)
! endif
-! if (itype(nres).eq.21) then
+! if (itype(nres,1).eq.21) then
! theta(nres)=90.0d0*deg2rad
! phi(nres)=180.0d0*deg2rad
! vbld(nres)=3.8d0
+(c(j,i+1)-c(j,i))*vbld_inv(i+1))
! in wham c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)
enddo
- iti=itype(i)
+ iti=itype(i,1)
di=dist(i,nres+i)
!#ifndef WHAM_RUN
! 10/03/12 Adam: Correction for zero SC-SC bond length
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1 .and. di.eq.0.0d0) &
- di=dsc(itype(i))
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
+ di=dsc(itype(i,1))
vbld(i+nres)=di
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
vbld_inv(i+nres)=1.0d0/di
else
vbld_inv(i+nres)=0.0d0
endif
if(me.eq.king.or..not.out1file)then
if (lprn) &
- write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),&
+ write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
rad2deg*alph(i),rad2deg*omeg(i)
endif
enddo
else if (lprn) then
do i=2,nres
- iti=itype(i)
+ iti=itype(i,1)
if(me.eq.king.or..not.out1file) &
- write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),&
+ write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,dist(i,i-1),&
rad2deg*theta(i),rad2deg*phi(i)
enddo
endif
enddo
enddo
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
do j=1,3
dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
enddo
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=itype(i)
+ it=itype(i,1)
if ((it.ne.10).and.(it.ne.ntyp1)) then
!el if (it.ne.10) then
enddo
if (lprn) then
do i=2,nres
- iti=itype(i)
+ iti=itype(i,1)
if(me.eq.king.or..not.out1file) &
- write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),&
+ write (iout,'(a3,i4,3f10.5)') restyp(iti,1),i,xxref(i),&
yyref(i),zzref(i)
enddo
endif
integer :: i,j,ires,nscat
real(kind=8),dimension(3,20) :: sccor
real(kind=8) :: sccmj
+! print *,"I am in sccenter",ires,nscat
do j=1,3
sccmj=0.0D0
do i=1,nscat
do i=1,nres-1
vbld(i+1)=vbl
vbld_inv(i+1)=1.0d0/vbld(i+1)
- vbld(i+1+nres)=dsc(itype(i+1))
- vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
+ vbld(i+1+nres)=dsc(itype(i+1,1))
+ vbld_inv(i+1+nres)=dsc_inv(itype(i+1,1))
! print *,vbld(i+1),vbld(i+1+nres)
enddo
return
do j=1,3
gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
+gloc(nres-2,icg)*dtheta(j,1,3)
- if(itype(2).ne.10) then
+ if(itype(2,1).ne.10) then
gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
gloc(ialph(2,1)+nside,icg)*domega(j,1,2)
endif
do j=1,3
gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
- if(itype(2).ne.10) then
+ if(itype(2,1).ne.10) then
gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
endif
- if(itype(3).ne.10) then
+ if(itype(3,1).ne.10) then
gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
endif
gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
dtheta(j,1,5)
- if(itype(3).ne.10) then
+ if(itype(3,1).ne.10) then
gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
endif
- if(itype(4).ne.10) then
+ if(itype(4,1).ne.10) then
gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
endif
+gloc(i-1,icg)*dphi(j,2,i+2)+ &
gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
gloc(nres+i-3,icg)*dtheta(j,1,i+2)
- if(itype(i).ne.10) then
+ if(itype(i,1).ne.10) then
gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
endif
- if(itype(i+1).ne.10) then
+ if(itype(i+1,1).ne.10) then
gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
+gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
endif
dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
+gloc(2*nres-6,icg)* &
dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
- if(itype(nres-2).ne.10) then
+ if(itype(nres-2,1).ne.10) then
gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
domega(j,2,nres-2)
endif
- if(itype(nres-1).ne.10) then
+ if(itype(nres-1,1).ne.10) then
gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
domega(j,1,nres-1)
do j=1,3
gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
gloc(2*nres-5,icg)*dtheta(j,2,nres)
- if(itype(nres-1).ne.10) then
+ if(itype(nres-1,1).ne.10) then
gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
domega(j,2,nres-1)
enddo
! The side-chain vector derivatives
do i=2,nres-1
- if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
+gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
! enddo
if (nres.lt.2) return
- if ((nres.lt.3).and.(itype(1).eq.10)) return
- if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+ if ((nres.lt.3).and.(itype(1,1).eq.10)) return
+ if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
do j=1,3
!c Derviative was calculated for oposite vector of side chain therefore
! there is "-" sign before gloc_sc
dtauangle(j,1,1,3)
gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
dtauangle(j,1,2,3)
- if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+ if ((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
gxcart(j,1)= gxcart(j,1) &
-gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
endif
enddo
endif
- if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1)) &
+ if ((nres.ge.3).and.(itype(3,1).ne.10).and.(itype(3,1).ne.ntyp1)) &
then
do j=1,3
gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
! Calculating the remainder of dE/ddc2
do j=1,3
- if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
- if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
+ if((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
+ if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
- if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1)) &
+ if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,1).ne.ntyp1)) &
then
gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
!c the - above is due to different vector direction
! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
endif
endif
- if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+ if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
! write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3)
endif
- if ((itype(3).ne.10).and.(nres.ge.3)) then
+ if ((itype(3,1).ne.10).and.(nres.ge.3)) then
gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
! write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
endif
- if ((itype(4).ne.10).and.(nres.ge.4)) then
+ if ((itype(4,1).ne.10).and.(nres.ge.4)) then
gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
! write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
endif
-! write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2"
+! write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
enddo
! If there are more than five residues
if(nres.ge.5) then
do i=3,nres-2
do j=1,3
! write(iout,*) "before", gcart(j,i)
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
*dtauangle(j,2,3,i+1) &
-gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
*dtauangle(j,1,2,i+2)
! write(iout,*) "new",j,i,
! & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
- if (itype(i-1).ne.10) then
+ if (itype(i-1,1).ne.10) then
gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
*dtauangle(j,3,3,i+1)
endif
- if (itype(i+1).ne.10) then
+ if (itype(i+1,1).ne.10) then
gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
*dtauangle(j,3,1,i+2)
gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
*dtauangle(j,3,2,i+2)
endif
endif
- if (itype(i-1).ne.10) then
+ if (itype(i-1,1).ne.10) then
gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
dtauangle(j,1,3,i+1)
endif
- if (itype(i+1).ne.10) then
+ if (itype(i+1,1).ne.10) then
gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
dtauangle(j,2,2,i+2)
! write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
! & dtauangle(j,2,2,i+2)
endif
- if (itype(i+2).ne.10) then
+ if (itype(i+2,1).ne.10) then
gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
dtauangle(j,2,1,i+3)
endif
! Setting dE/ddnres-1
if(nres.ge.4) then
do j=1,3
- if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
+ if ((itype(nres-1,1).ne.10).and.(itype(nres-1,1).ne.ntyp1)) then
gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
*dtauangle(j,2,3,nres)
! write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
! & dtauangle(j,2,3,nres), gxcart(j,nres-1)
- if (itype(nres-2).ne.10) then
+ if (itype(nres-2,1).ne.10) then
gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
*dtauangle(j,3,3,nres)
endif
- if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+ if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
*dtauangle(j,3,1,nres+1)
gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
*dtauangle(j,3,2,nres+1)
endif
endif
- if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
+ if ((itype(nres-2,1).ne.10).and.(itype(nres-2,1).ne.ntyp1)) then
gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
dtauangle(j,1,3,nres)
endif
- if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+ if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
dtauangle(j,2,2,nres+1)
! write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
-! & dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres)
+! & dtauangle(j,2,2,nres+1), itype(nres-1,1),itype(nres,1)
endif
enddo
endif
! Settind dE/ddnres
- if ((nres.ge.3).and.(itype(nres).ne.10).and. &
- (itype(nres).ne.ntyp1))then
+ if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
+ (itype(nres,1).ne.ntyp1))then
do j=1,3
gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
*dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
write (iout,100)
do i=1,nres
- write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+ write (iout,110) restyp(itype(i,1),1),i,c(1,i),c(2,i),&
c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
enddo
100 format (//' alpha-carbon coordinates ',&
DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
do i=nnt,nct
- if (itype(i).ne.10) then
-!d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
+ if (itype(i,1).ne.10) then
+!d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
fail=.true.
ii=0
do while (fail .and. ii .le. maxsi)
- call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
+ call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
ii = ii+1
enddo
endif
! include 'COMMON.BOUNDS'
! include 'COMMON.MD'
! include 'COMMON.SETUP'
- character(len=4),dimension(:),allocatable :: sequence !(maxres)
+ character(len=4),dimension(:,:),allocatable :: sequence !(maxres,maxmolecules)
! integer :: rescode
! double precision x(maxvar)
character(len=256) :: pdbfile
!-----------------------------
allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
- allocate(itype(maxres)) !(maxres)
+ allocate(itype(maxres,5)) !(maxres)
+ allocate(istype(maxres))
!
! Zero out tables.
!
c(:,:)=0.0D0
dc(:,:)=0.0D0
- itype(:)=0
+ itype(:,:)=0
!-----------------------------
!
! Body
!el if(.not.allocated(itype_pdb))
allocate(itype_pdb(nres))
do i=1,nres
- itype_pdb(i)=itype(i)
+ itype_pdb(i)=itype(i,1)
enddo
close (ipdbin)
nnt=nstart_sup
write(iout,*)'Adding sidechains'
maxsi=1000
do i=2,nres-1
- iti=itype(i)
- if (iti.ne.10 .and. itype(i).ne.ntyp1) then
+ iti=itype(i,1)
+ if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
nsi=0
fail=.true.
do while (fail.and.nsi.le.maxsi)
enddo
endif
endif
-
+
if (indpdb.eq.0) then
+ nres_molec(:)=0
+ allocate(sequence(maxres,5))
+
+ if (protein) then
! Read sequence if not taken from the pdb file.
- read (inp,*) nres
+ molec=1
+ read (inp,*) nres_molec(molec)
! print *,'nres=',nres
- allocate(sequence(nres))
if (iscode.gt.0) then
- read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
+ read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
else
- read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
+ read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
endif
! Convert sequence to numeric code
- do i=1,nres
- itype(i)=rescode(i,sequence(i),iscode)
+
+ do i=1,nres_molec(molec)
+ itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
+ enddo
+ endif
+ if (nucleic) then
+! Read sequence if not taken from the pdb file.
+ molec=2
+ read (inp,*) nres_molec(molec)
+! print *,'nres=',nres
+! allocate(sequence(maxres,5))
+! if (iscode.gt.0) then
+ read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
+! Convert sequence to numeric code
+
+ do i=1,nres_molec(molec)
+ istype(i)=sugarcode(sequence(i,molec)(1:1),i)
+ itype(i,1)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
enddo
+ endif
+
+ if (ions) then
+! Read sequence if not taken from the pdb file.
+ molec=5
+ read (inp,*) nres_molec(molec)
+! print *,'nres=',nres
+ read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
+! Convert sequence to numeric code
+ print *,nres_molec(molec)
+ do i=1,nres_molec(molec)
+ itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
+ enddo
+ endif
+ nres=0
+ do i=1,5
+ nres=nres+nres_molec(i)
+ print *,nres_molec(i)
+ enddo
+
! Assign initial virtual bond lengths
!elwrite(iout,*) "test_alloc"
if(.not.allocated(vbld)) allocate(vbld(2*nres))
vbld_inv(i)=vblinv
enddo
do i=2,nres-1
- vbld(i+nres)=dsc(iabs(itype(i)))
- vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
-! write (iout,*) "i",i," itype",itype(i),
-! & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
+ vbld(i+nres)=dsc(iabs(itype(i,1)))
+ vbld_inv(i+nres)=dsc_inv(iabs(itype(i,1)))
+! write (iout,*) "i",i," itype",itype(i,1),
+! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
enddo
endif
! print *,nres
-! print '(20i4)',(itype(i),i=1,nres)
+! print '(20i4)',(itype(i,1),i=1,nres)
!----------------------------
!el reallocate tables
! do i=1,maxres2
! enddo
! enddo
! do i=1,maxres
-!elwrite(iout,*) "itype",i,itype(i)
-! itype_alloc(i)=itype(i)
+!elwrite(iout,*) "itype",i,itype(i,1)
+! itype_alloc(i)=itype(i,1)
! enddo
! deallocate(c)
! enddo
! enddo
! do i=1,nres+2
-! itype(i)=itype_alloc(i)
+! itype(i,1)=itype_alloc(i)
! itel(i)=0
! enddo
!--------------------------
do i=1,nres
#ifdef PROCOR
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
#else
- if (itype(i).eq.ntyp1) then
+ if (itype(i,1).eq.ntyp1) then
#endif
itel(i)=0
#ifdef PROCOR
- else if (iabs(itype(i+1)).ne.20) then
+ else if (iabs(itype(i+1,1)).ne.20) then
#else
- else if (iabs(itype(i)).ne.20) then
+ else if (iabs(itype(i,1)).ne.20) then
#endif
itel(i)=1
else
if(me.eq.king.or..not.out1file)then
write (iout,*) "ITEL"
do i=1,nres-1
- write (iout,*) i,itype(i),itel(i)
+ write (iout,*) i,itype(i,1),itel(i)
enddo
print *,'Call Read_Bridge.'
endif
write (iout,'(a)') 'Boundaries in phi angle sampling:'
do i=1,nres
write (iout,'(a3,i5,2f10.1)') &
- restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
+ restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
enddo
#ifdef MP
endif
#endif
nct=nres
!d print *,'NNT=',NNT,' NCT=',NCT
- if (itype(1).eq.ntyp1) nnt=2
- if (itype(nres).eq.ntyp1) nct=nct-1
+ if (itype(1,1).eq.ntyp1) nnt=2
+ if (itype(nres,1).eq.ntyp1) nct=nct-1
if (pdbref) then
if(me.eq.king.or..not.out1file) &
write (iout,'(a,i3)') 'nsup=',nsup
nstart_seq=nnt
if (nsup.le.(nct-nnt+1)) then
do i=0,nct-nnt+1-nsup
- if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
+ if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
nstart_seq=nnt+i
goto 111
endif
stop
else
do i=0,nsup-(nct-nnt+1)
- if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) &
+ if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
then
nstart_sup=nstart_sup+i
nsup=nct-nnt+1
icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
enddo
if(me.eq.king.or..not.out1file) &
- write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',&
+ write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
icont_ref(1,i),' ',&
- restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
+ restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
enddo
endif
endif
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
dc(j,i+nres)=c(j,i+nres)-c(j,i)
dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
!elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
do i=2,nres-1
omeg(i)=-120d0*deg2rad
- if (itype(i).le.0) omeg(i)=-omeg(i)
+ if (itype(i,1).le.0) omeg(i)=-omeg(i)
enddo
else
if(me.eq.king.or..not.out1file) &
do i=1,nss
i1=ihpb(i)-nres
i2=jhpb(i)-nres
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,1)
+ it2=itype(i2,1)
if (me.eq.king.or..not.out1file) &
write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
- restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),&
+ restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
ebr,forcon(i)
enddo
write (iout,'(a)')
! Check whether the specified bridging residues are cystines.
do i=1,ns
write(iout,*) i,iss(i)
- if (itype(iss(i)).ne.1) then
+ if (itype(iss(i),1).ne.1) then
if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
'Do you REALLY think that the residue ',&
- restyp(itype(iss(i))),i,&
+ restyp(itype(iss(i),1),1),i,&
' can form a disulfide bridge?!!!'
write (*,'(2a,i3,a)') &
'Do you REALLY think that the residue ',&
- restyp(itype(iss(i))),i,&
+ restyp(itype(iss(i),1),1),i,&
' can form a disulfide bridge?!!!'
#ifdef MPI
call MPI_Finalize(MPI_COMM_WORLD,ierror)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
do j=1,3
dc(j,i+nres)=c(j,i+nres)-c(j,i)
dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
!model write (iunit,'(a5,i6)') 'MODEL',1
if (nhfrag.gt.0) then
do j=1,nhfrag
- iti=itype(hfrag(1,j))
- itj=itype(hfrag(2,j))
+ iti=itype(hfrag(1,j),1)
+ itj=itype(hfrag(2,j),1)
if (j.lt.10) then
write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
'HELIX',j,'H',j,&
- restyp(iti),hfrag(1,j)-1,&
- restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+ restyp(iti,1),hfrag(1,j)-1,&
+ restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
else
write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
'HELIX',j,'H',j,&
- restyp(iti),hfrag(1,j)-1,&
- restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+ restyp(iti,1),hfrag(1,j)-1,&
+ restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
endif
enddo
endif
do j=1,nbfrag
- iti=itype(bfrag(1,j))
- itj=itype(bfrag(2,j)-1)
+ iti=itype(bfrag(1,j),1)
+ itj=itype(bfrag(2,j)-1,1)
write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
'SHEET',1,'B',j,2,&
- restyp(iti),bfrag(1,j)-1,&
- restyp(itj),bfrag(2,j)-2,0
+ restyp(iti,1),bfrag(1,j)-1,&
+ restyp(itj,1),bfrag(2,j)-2,0
if (bfrag(3,j).gt.bfrag(4,j)) then
- itk=itype(bfrag(3,j))
- itl=itype(bfrag(4,j)+1)
+ itk=itype(bfrag(3,j),1)
+ itl=itype(bfrag(4,j)+1,1)
write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
'SHEET',2,'B',j,2,&
- restyp(itl),bfrag(4,j),&
- restyp(itk),bfrag(3,j)-1,-1,&
- "N",restyp(itk),bfrag(3,j)-1,&
- "O",restyp(iti),bfrag(1,j)-1
+ restyp(itl,1),bfrag(4,j),&
+ restyp(itk,1),bfrag(3,j)-1,-1,&
+ "N",restyp(itk,1),bfrag(3,j)-1,&
+ "O",restyp(iti,1),bfrag(1,j)-1
else
- itk=itype(bfrag(3,j))
- itl=itype(bfrag(4,j)-1)
+ itk=itype(bfrag(3,j),1)
+ itl=itype(bfrag(4,j)-1,1)
write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
'SHEET',2,'B',j,2,&
- restyp(itk),bfrag(3,j)-1,&
- restyp(itl),bfrag(4,j)-2,1,&
- "N",restyp(itk),bfrag(3,j)-1,&
- "O",restyp(iti),bfrag(1,j)-1
+ restyp(itk,1),bfrag(3,j)-1,&
+ restyp(itl,1),bfrag(4,j)-2,1,&
+ "N",restyp(itk,1),bfrag(3,j)-1,&
+ "O",restyp(iti,1),bfrag(1,j)-1
ichain=1
ires=0
do i=nnt,nct
- iti=itype(i)
- iti1=itype(i+1)
+ iti=itype(i,1)
+ iti1=itype(i+1,1)
if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
if (iti.eq.ntyp1) then
ichain=ichain+1
ires=ires+1
iatom=iatom+1
ica(i)=iatom
- write (iunit,10) iatom,restyp(iti),chainid(ichain),&
+ write (iunit,10) iatom,restyp(iti,1),chainid(ichain),&
ires,(c(j,i),j=1,3),vtot(i)
if (iti.ne.10) then
iatom=iatom+1
- write (iunit,20) iatom,restyp(iti),chainid(ichain),&
+ write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
ires,(c(j,nres+i),j=1,3),&
vtot(i+nres)
endif
enddo
write (iunit,'(a)') 'TER'
do i=nnt,nct-1
- if (itype(i).eq.ntyp1) cycle
- if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+ if (itype(i,1).eq.ntyp1) cycle
+ if (itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1) then
write (iunit,30) ica(i),ica(i+1)
- else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+ else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
write (iunit,30) ica(i),ica(i+1),ica(i)+1
- else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+ else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
write (iunit,30) ica(i),ica(i)+1
endif
enddo
- if (itype(nct).ne.10) then
+ if (itype(nct,1).ne.10) then
write (iunit,30) ica(nct),ica(nct)+1
endif
do i=1,nss
write (imol2,'(a)') '\@<TRIPOS>ATOM'
do i=nnt,nct
write (zahl,'(i3)') i
- pom=ucase(restyp(itype(i)))
+ pom=ucase(restyp(itype(i,1),1))
res_num = pom(:3)//zahl(2:)
write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
enddo
write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
do i=nnt,nct
write (zahl,'(i3)') i
- pom = ucase(restyp(itype(i)))
+ pom = ucase(restyp(itype(i,1),1))
res_num = pom(:3)//zahl(2:)
write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
enddo
write (iout,'(7a)') ' Res ',' d',' Theta',&
' Phi',' Dsc',' Alpha',' Omega'
do i=1,nres
- iti=itype(i)
- write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),&
+ iti=itype(i,1)
+ write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
rad2deg*omeg(i)
enddo
! write (iout,100)
! do i=1,nres
-! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
! enddo
! 100 format (//' alpha-carbon coordinates ',&
'inertia','Pstok'
write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
do i=1,ntyp
- write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
+ write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
do j=2,nbondterm(i)
write (iout,'(13x,3f10.5)') &
' ATHETA0 ',' A1 ',' A2 ',&
' B1 ',' B2 '
do i=1,ntyp
- write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
+ write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
enddo
write (iout,'(/a/9x,5a/79(1h-))') &
' ALPH0 ',' ALPH1 ',' ALPH2 ',&
' ALPH3 ',' SIGMA0C '
do i=1,ntyp
- write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
+ write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
(polthet(j,i),j=0,3),sigc0(i)
enddo
write (iout,'(/a/9x,5a/79(1h-))') &
' THETA0 ',' SIGMA0 ',' G1 ',&
' G2 ',' G3 '
do i=1,ntyp
- write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
+ write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
sig0(i),(gthet(j,i),j=1,3)
enddo
else
' theta0 ',' a1*10^2 ',' a2*10^2 ',&
' b1*10^1 ',' b2*10^1 '
do i=1,ntyp
- write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
+ write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
a0thet(i),(100*athet(j,i,1,1),j=1,2),&
(10*bthet(j,i,1,1),j=1,2)
enddo
' alpha0 ',' alph1 ',' alph2 ',&
' alhp3 ',' sigma0c '
do i=1,ntyp
- write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
+ write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
(polthet(j,i),j=0,3),sigc0(i)
enddo
write (iout,'(/a/9x,5a/79(1h-))') &
' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
' G2 ',' G3*10^1 '
do i=1,ntyp
- write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
+ write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
enddo
endif
nlobi=nlob(i)
if (nlobi.gt.0) then
if (LaTeX) then
- write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
+ write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
'log h',(bsc(j,i),j=1,nlobi)
call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
write (iout,'(/a)') 'One-body parameters:'
write (iout,'(a,4x,a)') 'residue','sigma'
- write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
+ write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
endif
! goto 50
!----------------------- LJK potential --------------------------------
call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
write (iout,'(/a)') 'One-body parameters:'
write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
- write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
+ write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
i=1,ntyp)
endif
! goto 50
write (iout,'(/a)') 'One-body parameters:'
write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
' chip ',' alph '
- write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
+ write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
chip(i),alp(i),i=1,ntyp)
endif
! goto 50
write (iout,'(/a)') 'One-body parameters:'
write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
's||/s_|_^2',' chip ',' alph '
- write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
+ write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
sigii(i),chip(i),alp(i),i=1,ntyp)
endif
case default
endif
if (lprint) then
write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
- restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
+ restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
endif
enddo
use control_data
use compare_data
use MPI_data
- use control, only: rescode
+ use control, only: rescode,sugarcode
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.LOCAL'
! include 'COMMON.CONTROL'
! include 'COMMON.DISTFIT'
! include 'COMMON.SETUP'
- integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
+ integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
! ishift_pdb
logical :: lprn=.true.,fail
real(kind=8),dimension(3) :: e1,e2,e3
character(len=5) :: atom
character(len=80) :: card
real(kind=8),dimension(3,20) :: sccor
- integer :: kkk,lll,icha,kupa !rescode,
+ integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
+ integer :: isugar
+ character*1 :: sugar
real(kind=8) :: cou
+ real(kind=8),dimension(3) :: ccc
!el local varables
integer,dimension(2,maxres/3) :: hfrag_alloc
integer,dimension(4,maxres/3) :: bfrag_alloc
real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
-
+ real(kind=8),dimension(:,:), allocatable :: c_temporary
+ integer,dimension(:,:) , allocatable :: itype_temporary
efree_temp=0.0d0
ibeg=1
ishift1=0
ishift=0
+ molecule=1
+ counter=0
! write (2,*) "UNRES_PDB",unres_pdb
ires=0
ires_old=0
else if (card(:3).eq.'TER') then
! End current chain
ires_old=ires+2
+ ishift=ishift+1
ishift1=ishift1+1
- itype(ires_old)=ntyp1
- itype(ires_old-1)=ntyp1
+ itype(ires_old,molecule)=ntyp1_molec(molecule)
+ itype(ires_old-1,molecule)=ntyp1_molec(molecule)
+ nres_molec(molecule)=nres_molec(molecule)+2
ibeg=2
! write (iout,*) "Chain ended",ires,ishift,ires_old
if (unres_pdb) then
! Read free energy
if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
! Fish out the ATOM cards.
+! write(iout,*) 'card',card(1:20)
if (index(card(1:4),'ATOM').gt.0) then
read (card(12:16),*) atom
! write (iout,*) "! ",atom," !",ires
ishift=ires-1
if (res.ne.'GLY' .and. res.ne. 'ACE') then
ishift=ishift-1
- itype(1)=ntyp1
+ itype(1,1)=ntyp1
+ nres_molec(molecule)=nres_molec(molecule)+1
endif
ires=ires-ishift+ishift1
ires_old=ires
ishift1=ishift1-1 !!!!!
! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
ires=ires-ishift+ishift1
+ print *,ires,ishift,ishift1
ires_old=ires
ibeg=0
else
ishift=ishift-(ires-ishift+ishift1-ires_old-1)
ires=ires-ishift+ishift1
ires_old=ires
- endif
+ endif
+! print *,'atom',ires,atom
if (res.eq.'ACE' .or. res.eq.'NHE') then
- itype(ires)=10
+ itype(ires,1)=10
+ else
+ if (atom.eq.'CA '.or.atom.eq.'N ') then
+ molecule=1
+ itype(ires,molecule)=rescode(ires,res,0,molecule)
+! nres_molec(molecule)=nres_molec(molecule)+1
else
- itype(ires)=rescode(ires,res,0)
+ molecule=2
+ itype(ires,molecule)=rescode(ires,res(3:4),0,molecule)
+! nres_molec(molecule)=nres_molec(molecule)+1
+ endif
endif
else
ires=ires-ishift+ishift1
if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+! print *,ires,ishift,ishift1
! write (iout,*) "backbone ",atom
#ifdef DEBUG
write (iout,'(2i3,2x,a,3f8.3)') &
- ires,itype(ires),res,(c(j,ires),j=1,3)
+ ires,itype(ires,1),res,(c(j,ires),j=1,3)
#endif
iii=iii+1
+ nres_molec(molecule)=nres_molec(molecule)+1
do j=1,3
sccor(j,iii)=c(j,ires)
enddo
-! write (*,*) card(23:27),ires,itype(ires)
+ else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
+ atom.eq."C2'" .or. atom.eq."C3'" &
+ .or. atom.eq."C4'" .or. atom.eq."O4'")) then
+ read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
+!c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
+ print *,ires,ishift,ishift1
+ counter=counter+1
+! iii=iii+1
+! do j=1,3
+! sccor(j,iii)=c(j,ires)
+! enddo
+ do j=1,3
+ c(j,ires)=c(j,ires)+ccc(j)/5.0
+ enddo
+ if (counter.eq.5) then
+! iii=iii+1
+ nres_molec(molecule)=nres_molec(molecule)+1
+! do j=1,3
+! sccor(j,iii)=c(j,ires)
+! enddo
+ counter=0
+ endif
+! print *, "ATOM",atom(1:3)
+ else if (atom(1:3).eq."C5'") then
+ read (card(19:19),'(a1)') sugar
+ isugar=sugarcode(sugar,ires)
+ if (ibeg.eq.1) then
+ istype(1)=isugar
+ else
+ istype(ires)=isugar
+ endif
+ if (unres_pdb) then
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+ else
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+! write (*,*) card(23:27),ires,itype(ires,1)
else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
atom.ne.'N' .and. atom.ne.'C' .and. &
atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
- atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+ atom.ne.'OXT' .and. atom(:2).ne.'3H' &
+ .and. atom.ne.'P '.and. &
+ atom(1:1).ne.'H' .and. &
+ atom.ne.'OP1' .and. atom.ne.'OP2 ') then
! write (iout,*) "sidechain ",atom
+! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
+ if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
+! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
iii=iii+1
read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
endif
endif
enddo
! Calculate dummy residue coordinates inside the "chain" of a multichain
! system
nres=ires
+ if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
+! print *,'I have', nres_molec(:)
+
+ do k=1,5
+ if (nres_molec(k).eq.0) cycle
do i=2,nres-1
-! write (iout,*) i,itype(i)
-! if (itype(i).eq.ntyp1) then
-! write (iout,*) "dummy",i,itype(i)
+! write (iout,*) i,itype(i,1)
+! if (itype(i,1).eq.ntyp1) then
+! write (iout,*) "dummy",i,itype(i,1)
! do j=1,3
! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
! c(j,i)=(c(j,i-1)+c(j,i+1))/2
! dc(j,i)=c(j,i)
! enddo
! endif
- if (itype(i).eq.ntyp1) then
- if (itype(i+1).eq.ntyp1) then
+ if (itype(i,k).eq.ntyp1_molec(k)) then
+ if (itype(i+1,k).eq.ntyp1_molec(k)) then
+ if (itype(i+2,k).eq.0) then
+ print *,"masz sieczke"
+ do j=1,5
+ if (itype(i+2,j).ne.0) then
+ itype(i+1,k)=0
+ itype(i+1,j)=ntyp1_molec(j)
+ nres_molec(k)=nres_molec(k)-1
+ nres_molec(j)=nres_molec(j)+1
+ go to 3331
+ endif
+ enddo
+ 3331 continue
+ endif
! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
+! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
if (unres_pdb) then
! 2/15/2013 by Adam: corrected insertion of the last dummy residue
! print *,i,'tu dochodze'
c(j,nres+i)=c(j,i)
enddo
endif !unres_pdb
- else !itype(i+1).eq.ntyp1
+ else !itype(i+1,1).eq.ntyp1
if (unres_pdb) then
! 2/15/2013 by Adam: corrected insertion of the first dummy residue
call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
c(j,nres+i)=c(j,i)
enddo
endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
+ endif !itype(i+1,1).eq.ntyp1
endif !itype.eq.ntyp1
enddo
+ enddo
! Calculate the CM of the last side chain.
if (iii.gt.0) then
if (unres_pdb) then
! nres=ires
nsup=nres
nstart_sup=1
- if (itype(nres).ne.10) then
+! print *,"molecule",molecule
+ if (itype(nres,1).ne.10) then
nres=nres+1
- itype(nres)=ntyp1
+ itype(nres,molecule)=ntyp1_molec(molecule)
+ nres_molec(molecule)=nres_molec(molecule)+1
if (unres_pdb) then
! 2/15/2013 by Adam: corrected insertion of the last dummy residue
call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
enddo
endif
endif
+! print *,'I have',nres, nres_molec(:)
+
!el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
#ifdef WHAM_RUN
if (nres.ne.nres0) then
c(j,nres+1)=c(j,1)
c(j,2*nres)=c(j,nres)
enddo
- if (itype(1).eq.ntyp1) then
+
+ if (itype(1,1).eq.ntyp1) then
nsup=nsup-1
nstart_sup=2
if (unres_pdb) then
enddo
endif
endif
+! First lets assign correct dummy to correct type of chain
+! 1) First residue
+ if (itype(1,1).eq.ntyp1) then
+ if (itype(2,1).eq.0) then
+ do j=2,5
+ if (itype(2,j).ne.0) then
+ itype(1,1)=0
+ itype(1,j)=ntyp1_molec(j)
+ nres_molec(1)=nres_molec(1)-1
+ nres_molec(j)=nres_molec(j)+1
+ go to 3231
+ endif
+ enddo
+3231 continue
+ endif
+ endif
+ print *,'I have',nres, nres_molec(:)
+
! Copy the coordinates to reference coordinates
! do i=1,2*nres
! do j=1,3
write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
"Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
do ires=1,nres
- write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
- restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
+ write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
+ (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
(c(j,ires+nres),j=1,3)
enddo
endif
"Backbone and SC coordinates as read from the PDB"
do ires=1,nres
write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
- ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
+ ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
(c(j,nres+ires),j=1,3)
enddo
endif
+! NOW LETS ROCK! SORTING
+ allocate(c_temporary(3,2*nres))
+ allocate(itype_temporary(nres,5))
+ itype_temporary(:,:)=0
+ seqalingbegin=1
+ do k=1,5
+ do i=1,nres
+ if (itype(i,k).ne.0) then
+ do j=1,3
+ c_temporary(j,seqalingbegin)=c(j,i)
+ c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
+
+ enddo
+ itype_temporary(seqalingbegin,k)=itype(i,k)
+ seqalingbegin=seqalingbegin+1
+ endif
+ enddo
+ enddo
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=c_temporary(j,i)
+ enddo
+ enddo
+ do k=1,5
+ do i=1,nres
+ itype(i,k)=itype_temporary(i,k)
+ enddo
+ enddo
+ if (itype(1,1).eq.ntyp1) then
+ nsup=nsup-1
+ nstart_sup=2
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,4)-c(j,3))/2.0
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ endif
+
+ if (lprn) then
+ write (iout,'(/a)') &
+ "Cartesian coordinates of the reference structure after sorting"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
+ "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
+ (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
+ (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+ print *,seqalingbegin,nres
if(.not.allocated(vbld)) then
allocate(vbld(2*nres))
do i=1,2*nres
lll=lll+1
!c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
if (i.gt.1) then
- if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
+ if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
chain_length=lll-1
kkk=kkk+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
do kkk=1,nperm
write (iout,*) "nowa struktura", nperm
do i=1,nres
- write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
+ write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
cref(2,i,kkk),&
cref(3,i,kkk),cref(1,nres+i,kkk),&
cref(2,nres+i,kkk),cref(3,nres+i,kkk)
write(iout,*) "shield_mode",shield_mode
!C Varibles set size of box
with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+ protein=index(controlcard,"PROTEIN").gt.0
+ ions=index(controlcard,"IONS").gt.0
+ nucleic=index(controlcard,"NUCLEIC").gt.0
write (iout,*) "with_theta_constr ",with_theta_constr
AFMlog=(index(controlcard,'AFM'))
selfguide=(index(controlcard,'SELFGUIDE'))
" stochastic forces of fully exposed sites"
write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
do i=1,ntyp
- write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
+ write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i),&
gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
enddo
endif
do ij=1,2
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
igall=igall+1
if (mask_side(i).eq.1) then
ig=ig+1
! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
! do i=1,nres
! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
-! restyp(itype(i)),i,(c(j,i),j=1,3),&
+! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
! (c(j,i+nres),j=1,3)
! enddo
!el----------------------------
enddo
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
IF (mask_side(i).eq.1) THEN
ig=ig+1
galphai=0.0D0
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
IF (mask_side(i).eq.1) THEN
ig=ig+1
gomegai=0.0D0
do ij=1,2
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
igall=igall+1
if (mask_side(i).eq.1) then
ig=ig+1
do ij=1,2
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
igall=igall+1
if (mask_side(i).eq.1) then
ig=ig+1
! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
! do i=1,nres
! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
-! restyp(itype(i)),i,(c(j,i),j=1,3),&
+! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
! (c(j,i+nres),j=1,3)
! enddo
!el----------------------------
sc_dist_cutoff=5.0D0
! Don't do glycine or ends
- i=itype(res_pick)
+ i=itype(res_pick,1)
if (i.eq.10 .or. i.eq.ntyp1) return
! Freeze everything (later will relax only selected side-chains)
!rc cur_e=orig_e
nres_moved=0
do i=2,nres-1
-! Don't do glycine (itype(j)==10)
- if (itype(i).ne.10) then
+! Don't do glycine (itype(j,1)==10)
+ if (itype(i,1).ne.10) then
sc_dist=dist(nres+i,nres+res_pick)
else
sc_dist=sc_dist_cutoff
n_try=0
do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
! Move the selected residue (don't worry if it fails)
- call gen_side(iabs(itype(res_pick)),theta(res_pick+1),&
+ call gen_side(iabs(itype(res_pick,1)),theta(res_pick+1),&
alph(res_pick),omeg(res_pick),fail)
! Minimize the side-chains starting from the new arrangement
! "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
! do i=1,nres
! write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
-! restyp(itype(i)),i,(c(j,i),j=1,3),&
+! restyp(itype(i,1)),i,(c(j,i),j=1,3),&
! (c(j,i+nres),j=1,3)
! enddo
! call etotal(energia)
enddo
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
IF (mask_side(i).eq.1) THEN
ig=ig+1
galphai=0.0D0
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
IF (mask_side(i).eq.1) THEN
ig=ig+1
gomegai=0.0D0
do ij=1,2
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
igall=igall+1
if (mask_side(i).eq.1) then
ig=ig+1
ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
- itypi1=iabs(itype(i+1))
+ itypi=iabs(itype(i,1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do j=istart(i,iint),iend(i,iint)
IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
dscj_inv=dsc_inv(itypj)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
projects / unres4.git / commitdiff