!
allocate(dsc(ntyp1)) !(ntyp1)
allocate(dsc_inv(ntyp1)) !(ntyp1)
+ allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
+ allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+ allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
allocate(nbondterm(ntyp)) !(ntyp)
allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
- allocate(msc(ntyp+1)) !(ntyp+1)
- allocate(isc(ntyp+1)) !(ntyp+1)
- allocate(restok(ntyp+1)) !(ntyp+1)
+ allocate(msc(ntyp+1,5)) !(ntyp+1)
+ allocate(isc(ntyp+1,5)) !(ntyp+1)
+ allocate(restok(ntyp+1,5)) !(ntyp+1)
allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
allocate(long_r_sidechain(ntyp))
allocate(short_r_sidechain(ntyp))
endif
enddo
#else
- read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
- do i=1,ntyp
+ read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
+ do i=1,ntyp_molec(1)
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
- j=1,nbondterm(i)),msc(i),isc(i),restok(i)
+ j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
dsc(i) = vbldsc0(1,i)
if (i.eq.10) then
dsc_inv(i)=0.0D0
endif
enddo
#endif
+ read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
+ do i=1,ntyp_molec(2)
+ nbondterm_nucl(i)=1
+ read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
+! dsc(i) = vbldsc0_nucl(1,i)
+! if (i.eq.10) then
+! dsc_inv(i)=0.0D0
+! else
+! dsc_inv(i)=1.0D0/dsc(i)
+! endif
+ enddo
+! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
+! do i=1,ntyp_molec(2)
+! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
+! aksc_nucl(j,i),abond0_nucl(j,i),&
+! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
+! dsc(i) = vbldsc0(1,i)
+! if (i.eq.10) then
+! dsc_inv(i)=0.0D0
+! else
+! dsc_inv(i)=1.0D0/dsc(i)
+! endif
+! enddo
+
if (lprint) then
write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
'inertia','Pstok'
- write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
+ write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
do i=1,ntyp
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)
+ vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
do j=2,nbondterm(i)
write (iout,'(13x,3f10.5)') &
vbldsc0(j,i),aksc(j,i),abond0(j,i)
close (ithep_pdb)
#endif
close(ithep)
+!--------------- Reading theta parameters for nucleic acid-------
+ read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
+ ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
+ nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
+ allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ allocate(aa0thet_nucl(nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxtheterm,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+
+ read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
+
+ aa0thet_nucl(:,:,:)=0.0d0
+ aathet_nucl(:,:,:,:)=0.0d0
+ bbthet_nucl(:,:,:,:,:)=0.0d0
+ ccthet_nucl(:,:,:,:,:)=0.0d0
+ ddthet_nucl(:,:,:,:,:)=0.0d0
+ eethet_nucl(:,:,:,:,:)=0.0d0
+ ffthet_nucl(:,:,:,:,:,:)=0.0d0
+ ggthet_nucl(:,:,:,:,:,:)=0.0d0
+
+ do i=1,nthetyp_nucl
+ do j=1,nthetyp_nucl
+ do k=1,nthetyp_nucl
+ read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
+ read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
+ read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
+ read (ithep_nucl,*,end=111,err=111) &
+ (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
+ read (ithep_nucl,*,end=111,err=111) &
+ (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
+ ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
+ llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
+ enddo
+ enddo
+ enddo
!-------------------------------------------
allocate(nlob(ntyp1)) !(ntyp1)
enddo
endif
#endif
+ allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
+ print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
+!el from energy module---------
+ allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+ allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+ allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
+ allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
+ allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
+ allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el---------------------------
+ nterm_nucl(:,:)=0
+ nlor_nucl(:,:)=0
+!el--------------------
+ read (itorp_nucl,*,end=113,err=113) &
+ (itortyp_nucl(i),i=1,ntyp_molec(2))
+ print *,itortyp_nucl(:)
+!c write (iout,*) 'ntortyp',ntortyp
+ do i=1,ntortyp_nucl
+ do j=1,ntortyp_nucl
+ read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
+ print *,nterm_nucl(i,j),nlor_nucl(i,j)
+ v0ij=0.0d0
+ si=-1.0d0
+ do k=1,nterm_nucl(i,j)
+ read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
+ v0ij=v0ij+si*v1_nucl(k,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_nucl(i,j)
+ read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
+ vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
+ v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
+ enddo
+ v0_nucl(i,j)=v0ij
+ enddo
+ enddo
+
! Read of Side-chain backbone correlation parameters
! Modified 11 May 2012 by Adasko
!CC
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,molecprev
+ 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,1)=ntyp1
- itype(ires_old-1,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
if (res.ne.'GLY' .and. res.ne. 'ACE') then
ishift=ishift-1
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,1)=10
else
- itype(ires,1)=rescode(ires,res,0,1)
+ 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
+ molecule=2
+ itype(ires,molecule)=rescode(ires,res(2: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,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
+ 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 '.and. atom.ne.'OP3'&
+ ) 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
+ else if ((ions).and.(card(1:6).eq.'HETATM')) then
+
+ read (card(12:16),*) atom
+ print *,"HETATOM", atom
+ read (card(18:20),'(a3)') res
+ if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
+ (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
+ .or.(atom(1:2).eq.'K ')) &
+ then
+ ires=ires+1
+ if (molecule.ne.5) molecprev=molecule
+ molecule=5
+ nres_molec(molecule)=nres_molec(molecule)+1
+ itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+ endif
+ endif !atom
enddo
10 write (iout,'(a,i5)') ' Number of residues found: ',ires
if (ires.eq.0) return
! 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,4 ! ions are without dummy
+ if (nres_molec(k).eq.0) cycle
do i=2,nres-1
! write (iout,*) i,itype(i,1)
! if (itype(i,1).eq.ntyp1) then
! dc(j,i)=c(j,i)
! enddo
! endif
- if (itype(i,1).eq.ntyp1) then
- if (itype(i+1,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,1).eq.ntyp1)=true
! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
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,1).ne.10) then
+! print *,"molecule",molecule
+ if ((itype(nres,1).ne.10)) then
nres=nres+1
- itype(nres,1)=ntyp1
+ if (molecule.eq.5) molecule=molecprev
+ 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,1).eq.ntyp1) then
nsup=nsup-1
nstart_sup=2
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,1),1),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
(c(j,nres+ires),j=1,3)
enddo
endif
+! NOW LETS ROCK! SORTING
+ allocate(c_temporary(3,2*nres))
+ allocate(itype_temporary(nres,5))
+ allocate(molnum(nres))
+ 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)
+ molnum(seqalingbegin)=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
if(me.eq.king.or..not.out1file) &
write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
eta
- gamp=scal_fric*(pstok+rwat)*eta
+ gamp=scal_fric*(pstok(i)+rwat)*eta
stdfp=dsqrt(2*Rb*t_bath/d_time)
allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
do i=1,ntyp
- gamsc(i)=scal_fric*(restok(i)+rwat)*eta
+ gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta
stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
enddo
if(me.eq.king.or..not.out1file)then
" 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,1),restok(i),&
+ write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
enddo
endif
! Get parameter filenames and open the parameter files.
call getenv_loc('BONDPAR',bondname)
open (ibond,file=bondname,status='old',readonly,shared)
+ call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old',readonly,shared)
call getenv_loc('ROTPAR',rotname)
open (ielep,file=elename,status='old',readonly,shared)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly,shared)
+
+ call getenv_loc('THETPAR_NUCL',thetname_nucl)
+ open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
+ call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+ open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
+ call getenv_loc('TORPAR_NUCL',torname_nucl)
+ open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
+ call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+ open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
+ call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+ open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
+
+
#elif (defined CRAY) || (defined AIX)
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
action='read')
! Get parameter filenames and open the parameter files.
call getenv_loc('BONDPAR',bondname)
open (ibond,file=bondname,status='old',action='read')
+ call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old',action='read')
+
! print *,"Processor",myrank," opened file IBOND"
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old',action='read')
! print *,"Processor",myrank," opened file IELEP"
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',action='read')
+
+ call getenv_loc('THETPAR_NUCL',thetname_nucl)
+ open (ithep_nucl,file=thetname_nucl,status='old',action='read')
+ call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+ open (irotam_nucl,file=rotname_nucl,status='old',action='read')
+ call getenv_loc('TORPAR_NUCL',torname_nucl)
+ open (itorp_nucl,file=torname_nucl,status='old',action='read')
+ call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+ open (itordp_nucl,file=tordname_nucl,status='old',action='read')
+ call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+ open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+
call getenv_loc('LIPTRANPAR',liptranname)
open (iliptranpar,file=liptranname,status='old',action='read')
call getenv_loc('TUBEPAR',tubename)
! Get parameter filenames and open the parameter files.
call getenv_loc('BONDPAR',bondname)
open (ibond,file=bondname,status='old')
+ call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old')
+
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old')
call getenv_loc('ROTPAR',rotname)
open (ielep,file=elename,status='old')
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old')
+
+ open (ithep_nucl,file=thetname_nucl,status='old')
+ call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+ open (irotam_nucl,file=rotname_nucl,status='old')
+ call getenv_loc('TORPAR_NUCL',torname_nucl)
+ open (itorp_nucl,file=torname_nucl,status='old')
+ call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+ open (itordp_nucl,file=tordname_nucl,status='old')
+ call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+ open (isidep_nucl,file=sidename_nucl,status='old')
+
call getenv_loc('LIPTRANPAR',liptranname)
open (iliptranpar,file=liptranname,status='old')
call getenv_loc('TUBEPAR',tubename)
! Get parameter filenames and open the parameter files.
call getenv_loc('BONDPAR',bondname)
open (ibond,file=bondname,status='old',action='read')
+ call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old',action='read')
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old',action='read')
call getenv_loc('ROTPAR',rotname)
open (ielep,file=elename,status='old',readonly)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly)
+
+ call getenv_loc('THETPAR_NUCL',thetname_nucl)
+ open (ithep_nucl,file=thetname_nucl,status='old',action='read')
+ call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+ open (irotam_nucl,file=rotname_nucl,status='old',action='read')
+ call getenv_loc('TORPAR_NUCL',torname_nucl)
+ open (itorp_nucl,file=torname_nucl,status='old',action='read')
+ call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+ open (itordp_nucl,file=tordname_nucl,status='old',action='read')
+ call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+ open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+
call getenv_loc('LIPTRANPAR',liptranname)
open (iliptranpar,file=liptranname,status='old',action='read')
call getenv_loc('TUBEPAR',tubename)