if (max_cache_traj_use.ne.1) &
print *,itime,"processor ",me," over cache ",ntwx_cache
do i=1,ntwx_cache-1
-
+ call returnbox
totT_cache(i)=totT_cache(i+1)
EK_cache(i)=EK_cache(i+1)
potE_cache(i)=potE_cache(i+1)
ugamma_cache(i,ntwx_cache)=ugamma(i)
uscdiff_cache(i,ntwx_cache)=uscdiff(i)
enddo
- call returnbox
+! call returnbox
do i=1,nres*2
do j=1,3
c_cache(j,i,ntwx_cache)=c(j,i)
nside=0
do i=2,nres-1
mnum=molnum(i)
+ write(iout,*) "i",molnum(i)
#ifdef WHAM_RUN
if (itype(i,1).ne.10) then
#else
"WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",&
"WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",&
"WCATPROT ","WCATCAT ","WNULL ","WNULL ","WNULL ",&
- "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPH O "/)
+ "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "/)
integer :: nprint_ene = 21
integer,dimension(n_ene) :: print_order = &
call esb(esbloc)
call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
else
+ etors_nucl=0.0d0
estr_nucl=0.0d0
ebe_nucl=0.0d0
evdwsb=0.0d0
!c
!c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!c
- print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+! print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
do i=iatel_s_nucl,iatel_e_nucl
if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
dxi=dc(1,i)
evdwij=evdwij*eps2rt
evdwsb=evdwsb+evdwij
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi,2),i,restyp(itypj,2),j, &
epsi,sigm,chi1,chi2,chip1,chip2, &
r012 = r06**2
k0 = 332.0*(2.0*2.0)/80.0
itmp=0
+
do i=1,4
itmp=itmp+nres_molec(i)
enddo
+ write(iout,*) "itmp",itmp
do i=itmp+1,itmp+nres_molec(5)-1
xi=c(1,i)
yi=c(2,i)
zi=c(3,i)
+
xi=mod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
yi=mod(yi,boxysize)
if (yj.lt.0) yj=yj+boxysize
zj=dmod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ write(iout,*) c(1,i),xi,xj,"xy",boxxsize
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
enddo
+ write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
enddo
enddo
allareout=1
!C change suggested by Ana -end
do i=1,nres-1
+ mnum=molnum(i)
if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
.and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
chain_end=i
c(1,i)=dmod(c(1,i),boxxsize)
c(2,i)=dmod(c(2,i),boxysize)
c(3,i)=dmod(c(3,i),boxzsize)
+ c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
+ c(2,i+nres)=dmod(c(2,i+nres),boxysize)
+ c(3,i+nres)=dmod(c(3,i+nres),boxzsize)
endif
enddo
return
((c(l,k+nres),l=1,3),k=nnt,nct)
write (iout,*) "Exit READ_CART"
write (iout,'(8f10.5)') &
- ((c(l,k),l=1,3),k=1,nres),&
+ ((c(l,k),l=1,3),k=1,nres)
+ write (iout,'(8f10.5)') &
((c(l,k+nres),l=1,3),k=nnt,nct)
call int_from_cart1(.true.)
write (iout,*) "Finish INT_TO_CART"
print *,"ENTER READ"
! Read bridging residues.
read (inp,*) ns
- print *,"ns",ns
+ write(iout,*) "ns",ns
+ call flush(iout)
if (ns.gt.0) then
allocate(iss(ns))
read (inp,*) (iss(i),i=1,ns)
! write (2,*) "UNRES_PDB",unres_pdb
ires=0
ires_old=0
+#ifdef WHAM_RUN
+ do i=1,nres
+ do j=1,5
+ itype(i,j)=0
+ enddo
+ enddo
+#endif
nres=0
iii=0
lsecondary=.false.
nhfrag=0
nbfrag=0
+ do j=1,5
+ nres_molec(j)=0
+ enddo
+
+
!-----------------------------
allocate(hfrag(2,maxres/3)) !(2,maxres/3)
allocate(bfrag(4,maxres/3)) !(4,maxres/3)
-
+ if(.not. allocated(istype)) allocate(istype(maxres))
do i=1,100000
read (ipdbin,'(a80)',end=10) card
! write (iout,'(a)') card
! nres_molec(molecule)=nres_molec(molecule)+1
else
molecule=2
- itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
+ itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
! nres_molec(molecule)=nres_molec(molecule)+1
read (card(19:19),'(a1)') sugar
isugar=sugarcode(sugar,ires)
molecule=5
nres_molec(molecule)=nres_molec(molecule)+1
print *,"HERE",nres_molec(molecule)
- itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
+ res=res(2:3)//' '
+ itype(ires,molecule)=rescode(ires,res,0,molecule)
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
endif
endif !atom
! NOW LETS ROCK! SORTING
allocate(c_temporary(3,2*nres))
allocate(itype_temporary(nres,5))
- allocate(molnum(nres+1))
+ if (.not.allocated(molnum)) allocate(molnum(nres+1))
+ if (.not.allocated(istype)) write(iout,*) &
+ "SOMETHING WRONG WITH ISTYTPE"
allocate(istype_temp(nres))
itype_temporary(:,:)=0
seqalingbegin=1
integer :: j,k
call alloc_compare_arrays
- if (indpdb.eq.0) then
+ if ((indpdb.eq.0).and.(.not.read_cart)) then
call chainbuild
write(iout,*) 'Warning: Calling chainbuild'
endif
read (intin,'(i5)',end=1100,err=1100) iconf
call read_angles(intin,*11)
call geom_to_var(nvar,varia)
- write(iout,*) 'Warning: Calling chainbuild'
+ write(iout,*) 'Warning: Calling chainbuild1'
call chainbuild
endif
write (iout,'(a,i7)') 'Conformation #',iconf
read (intin,'(i5)',end=11,err=11) iconf
call read_angles(intin,*11)
call geom_to_var(nvar,varia)
- write(iout,*) 'Warning: Calling chainbuild'
+ write(iout,*) 'Warning: Calling chainbuild2'
call chainbuild
endif
write (iout,'(a,i7)') 'Conformation #',iconf
! print *,'result received from worker ',man,' sending now'
call var_to_geom(nvar,varia)
- write(iout,*) 'Warning: Calling chainbuild'
+ write(iout,*) 'Warning: Calling chainbuild3'
call chainbuild
call etotal(energy_)
iconf=ind(2)
read (intin,'(i5)',end=1101,err=1101) iconf
call read_angles(intin,*11)
call geom_to_var(nvar,varia)
- write(iout,*) 'Warning: Calling chainbuild'
+ write(iout,*) 'Warning: Calling chainbuild4'
call chainbuild
endif
n=n+1
CG_COMM,muster,ierr)
call var_to_geom(nvar,varia)
- write(iout,*) 'Warning: Calling chainbuild'
+ write(iout,*) 'Warning: Calling chainbuild5'
call chainbuild
call etotal(energy_)
iconf=ind(2)
read (intin,'(i5)',end=11,err=11) iconf
call read_angles(intin,*11)
call geom_to_var(nvar,varia)
- write(iout,*) 'Warning: Calling chainbuild'
+ write(iout,*) 'Warning: Calling chainbuild5'
call chainbuild
endif
write (iout,'(a,i7)') 'Conformation #',iconf
# Set compiler flags for different sourcefiles
#================================================
if (Fortran_COMPILER_NAME STREQUAL "ifort")
- set (CMAKE_Fortran_FLAGS_RELEASE " ")
+ set (CMAKE_Fortran_FLAGS_RELEASE " -CB -g")
set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g ")
- set(FFLAGS0 "-fpp -mcmodel=medium -shared-intel " )
+ set(FFLAGS0 "-fpp -mcmodel=medium -shared-intel " )
elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
set(FFLAGS0 "-fpp -std=legacy -mcmodel=medium -g ")
elseif (Fortran_COMPILER_NAME STREQUAL "pgf90")
use calc_data
use geometry_data, only:c,dc,dc_norm
- use energy_data, only:itype,maxcont
+ use energy_data, only:itype,maxcont,molnum
! implicit none
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
! include 'COMMON.CALC'
! include 'COMMON.CONTPAR'
! include 'COMMON.LOCAL'
- integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2
+ integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2,mnum,mnum2
real(kind=8) :: csc !el,dist
real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,&
ddsc,ddla,ddlb
- integer :: ncont
+ integer :: ncont,mhum
integer,dimension(2,maxcont) :: icont
real(kind=8) :: u,v,a(3),b(3),dla,dlb
logical :: lprint
if (lprint) then
do i=1,nres
mnum=molnum(i)
- write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+ write (iout,110) restyp(itype(i,mnum),mnum),i,c(1,i),c(2,i),&
c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),&
dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
enddo
endif
110 format (a,'(',i3,')',9f8.3)
do i=ist,ien-kkk
- iti=iabs(itype(i))
- if (iti.le.0 .or. iti.gt.ntyp) cycle
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (iti.le.0 .or. iti.gt.ntyp_molec(mnum)) cycle
do j=i+kkk,ien
- itj=iabs(itype(j))
- if (itj.le.0 .or. itj.gt.ntyp) cycle
+ mnum2=molnum(j)
+ itj=iabs(itype(j,mnum2))
+ if (itj.le.0 .or. itj.gt.ntyp_molec(mnum2)) cycle
itypi=iti
itypj=itj
xj = c(1,nres+j)-c(1,nres+i)
if (lprint) then
write (iout,'(a)') 'Contact map:'
do i=1,ncont
+ mnum=molnum(i)
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,molnum(i1))
+ it2=itype(i2,molnum(i2))
+! print *,"CONTACT",i1,mnum,it1,it2
write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') &
- i,restyp(it1),i1,restyp(it2),i2,cscore(i),&
+ i,restyp(it1,mnum),i1,restyp(it2,mnum),i2,cscore(i),&
sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),&
omt1(i),omt2(i),omt12(i)
enddo
kkk=3
! print *,'nnt=',nnt,' nct=',nct
do i=nnt+kkk,nct
- iti=iabs(itype(i))
+ mnum=molnum(i)
+ iti=iabs(itype(i,1))
do j=nnt,i-kkk
- itj=iabs(itype(j))
+ mnum2=molnum(j)
+ itj=iabs(itype(j,1))
if (ipot.ne.4) then
! rcomp=sigmaii(iti,itj)+1.0D0
rcomp=facont*sigmaii(iti,itj)
if (lprint) then
write (iout,'(a)') 'Contact map:'
do i=1,ncont
+ mnum=molnum(i)
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,mnum),i1,restyp(it2,mnum),i2
enddo
endif
return
subroutine pept_cont(lprint,ncont,icont)
use geometry_data, only:c
- use energy_data, only:maxcont,nnt,nct,itype
+ use energy_data, only:maxcont,nnt,nct,itype,molnum
! implicit none
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
! include 'COMMON.FFIELD'
! include 'COMMON.NAMES'
integer :: ncont,icont(2,maxcont)
- integer :: i,j,k,kkk,i1,i2,it1,it2
+ integer :: i,j,k,kkk,i1,i2,it1,it2,mnum
logical :: lprint
!el real(kind=8) :: dist
real(kind=8) :: rcomp=5.5d0
if (lprint) then
write (iout,'(a)') 'PP contact map:'
do i=1,ncont
+ mnum=molnum(i)
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
enddo
endif
return
subroutine contacts_between_fragments(lprint,is,ncont,icont,&
ncont_interfrag,icont_interfrag)
- use energy_data, only:itype,maxcont
+ use energy_data, only:itype,maxcont,molnum
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
! include 'DIMENSIONS.COMPAR'
integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),&
icont_interfrag(2,maxcont,mmaxfrag)
logical :: OK1,OK2,lprint
- integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2
+ integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2,mnum
! Determine the contacts that occur within a fragment and between fragments.
do i=1,nfrag(1)
do j=1,i
do k=1,ncont_interfrag(ind)
i1=icont_interfrag(1,k,ind)
i2=icont_interfrag(2,k,ind)
- it1=itype(i1)
- it2=itype(i2)
+ mnum=molnum(i1)
+ it1=itype(i1,mnum)
+ it2=itype(i2,molnum(i2))
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
enddo
enddo
write (iout,*)
do k=1,ncont_interfrag(ind)
i1=icont_interfrag(1,k,ind)
i2=icont_interfrag(2,k,ind)
- it1=itype(i1)
- it2=itype(i2)
+ mnum=molnum(i1)
+ it1=itype(i1,mnum)
+ it2=itype(i2,molnum(i2))
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
enddo
enddo
enddo
subroutine elecont(lprint,ncont,icont,ist,ien)
use geometry_data, only:c
- use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2
+ use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2,molnum
! implicit none
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc
real(kind=8) :: elcutoff=-0.3d0
real(kind=8) :: elecutoff_14=-0.5d0
- integer :: ncont,icont(2,maxcont)
+ integer :: ncont,icont(2,maxcont),mnum
real(kind=8) :: econt(maxcont)
!
! Load the constants of peptide bond - peptide bond interactions.
do i=1,ncont
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,molnum(i1))
+ it2=itype(i2,molnum(i1))
write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
- i,restyp(it1),i1,restyp(it2),i2,econt(i)
+ i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i1)),i2,econt(i)
enddo
endif
! For given residues keep only the contacts with the greatest energy.
write (iout,*)
write (iout,*) 'Electrostatic contacts after pruning: '
do i=1,ncont
+ mnum=molnum(i)
i1=icont(1,i)
i2=icont(2,i)
- it1=itype(i1)
- it2=itype(i2)
+ mnum=molnum(i1)
+ it1=itype(i1,mnum)
+ it2=itype(i2,molnum(i2))
write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
- i,restyp(it1),i1,restyp(it2),i2,econt(i)
+ i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2,econt(i)
enddo
endif
return
use geometry_data, only:cref,nperm
use control_data, only:symetr
- use energy_data, only:nnt,nct,itype
+ use energy_data, only:nnt,nct,itype,molnum
! implicit none
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
(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,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
qq = qq+qqij+qqijCM
if (lprn) then
write (iout,*) "il",il," jl",jl,&
- " itype",itype(il),itype(jl)
+ " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
" dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
endif
(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,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
qq = qq+qqij+qqijCM
if (lprn) then
write (iout,*) "i",i," j",j," il",il," jl",jl,&
- " itype",itype(il),itype(jl)
+ " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
" dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
endif
(cref(3,kl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,kl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
- if (itype(il).ne.10 .or. itype(kl).ne.10) then
+ if (itype(il,molnum(il)).ne.10 .or. itype(kl,molnum(kl)).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
qq = qq+qqij+qqijCM
if (lprn) then
write (iout,*) "i",i," j",j," k",k," l",l," il",il,&
- " kl",kl," itype",itype(il),itype(kl)
+ " kl",kl," itype",itype(il,molnum(il)), &
+ itype(kl,molnum(kl))
write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",&
d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
endif
use control_data, only:symetr
use geometry_data, only:nperm,cref,c
- use energy_data, only:itype
+ use energy_data, only:itype,molnum
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
do kkk=1,nperm
nnsup=0
do i=1,nres
- if (itype(i).ne.ntyp1) then
+ if (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then
nnsup=nnsup+1
do j=1,3
cc(j,nnsup)=c(j,i)
subroutine define_fragments
use geometry_data, only:rad2deg
- use energy_data, only:itype
+ use energy_data, only:itype,molnum
use compare_data, only:nhfrag,nbfrag,bfrag,hfrag
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.INTERACT'
! include 'COMMON.NAMES'
integer :: nstrand,istrand(2,nres/2)
- integer :: nhairp,ihairp(2,nres/5)
+ integer :: nhairp,ihairp(2,nres/5),mnum
character(len=16) :: strstr(4)=reshape((/'helix','hairpin',&
'strand','strand pair'/),shape(strstr))
integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4
write (iout,*) "The following primary fragments were found:"
write (iout,*) "Helices:",nhfrag
do i=1,nhfrag
+ mnum=molnum(i)
i1=ifrag(1,1,i)
i2=ifrag(2,1,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
enddo
write (iout,*) "Hairpins:",nhairp
do i=nhfrag+1,nhfrag+nhairp
i1=ifrag(1,1,i)
i2=ifrag(2,1,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,molnum(i1))
+ it2=itype(i2,molnum(i2))
write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2
enddo
write (iout,*) "Far strand pairs:",nbfrag
do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
i1=ifrag(1,1,i)
i2=ifrag(2,1,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,molnum(i1))
+ it2=itype(i2,molnum(i1))
i3=ifrag(1,2,i)
i4=ifrag(2,2,i)
- it3=itype(i3)
- it4=itype(i4)
+ it3=itype(i3,molnum(i3))
+ it4=itype(i4,molnum(i4))
write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2,&
- restyp(it3),i3,restyp(it4),i4
+ i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2,&
+ restyp(it3,molnum(i3)),i3,restyp(it4,molnum(i4)),i4
enddo
write (iout,*) "Strands:",nstrand
do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
+ mnum=molnum(i)
i1=ifrag(1,1,i)
i2=ifrag(2,1,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(i3,2x,a,i4,2x,a,i4)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
enddo
call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),&
istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),&
nc_fragm(1,1),nc_req_setf(1,1))
write (iout,*) "Fragments after sorting:"
do i=1,nfrag(1)
+ mnum=molnum(i)
i1=ifrag(1,1,i)
i2=ifrag(2,1,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(i3,2x,a,i4,2x,a,i4,$)') &
- i,restyp(it1),i1,restyp(it2),i2
+ i,restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2
if (npiece(i,1).eq.1) then
write (iout,'(2x,a)') strstr(istruct(i))
else
i1=ifrag(1,2,i)
i2=ifrag(2,2,i)
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(2x,a,i4,2x,a,i4,2x,a)') &
- restyp(it1),i1,restyp(it2),i2,strstr(istruct(i))
+ restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2,strstr(istruct(i))
endif
enddo
return
subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup
- use energy_data, only:itype,maxcont
+ use energy_data, only:itype,maxcont,molnum
use compare_data, only:bfrag,hfrag,nbfrag,nhfrag
use compare, only:freeres
! implicit real*8 (a-h,o-z)
isecstr(nres)
logical :: lprint,lprint_sec,not_done !el,freeres
integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist
- integer :: nstrand,nbeta,nhelix,iii1,jjj1
+ integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum
real(kind=8) :: p1,p2
!rel external freeres
character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec))
write (iout,*)
write (iout,*) "Secondary structure"
do i=1,nres,80
+ mnum=molnum(i)
ist=i
ien=min0(i+79,nres)
write (iout,*)
write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
- write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien)
+ write (iout,'(80a1)') (onelet(itype(k,mnum)),k=ist,ien)
write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien)
enddo
write (iout,*)
write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
! call enerprint(energia(0),fT)
call enerprint(energia(0))
- write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
+ write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49)
write (iout,*) "ftors",ftors
!el call briefout(i,energia(0))
temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
endif
endif
potE(iii+1,iparm)=energia(0)
- do k=1,21
+ do k=1,49
enetb(k,iii+1,iparm)=energia(k)
enddo
! write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
#ifdef MPI
integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
#endif
- integer :: j,k,l,ii,itj,iprint
+ integer :: j,k,l,ii,itj,iprint,mnum,mnum1
if (.not.check_conf) then
conf_check=.true.
return
endif
call int_from_cart1(.false.)
do j=nnt+1,nct
- if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. &
+ mnum=molnum(j)
+ mnum1=molnum(j-1)
+ if (mnum.ne.1) cycle
+ if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
(vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
if (iprint.gt.0) &
write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
endif
enddo
do j=nnt,nct
- itj=itype(j)
- if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. &
+ mnum=molnum(j)
+ if (mnum.gt.3) cycle
+ itj=itype(j,mnum)
+ if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. &
(vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
if (iprint.gt.0) &
write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
! Store sidechain parameters
do i=1,ntyp
do j=1,ntyp
- aa_all(j,i,iparm)=aa(j,i)
- bb_all(j,i,iparm)=bb(j,i)
+ aa_all(j,i,iparm)=aa_aq(j,i)
+ bb_all(j,i,iparm)=bb_aq(j,i)
r0_all(j,i,iparm)=r0(j,i)
sigma_all(j,i,iparm)=sigma(j,i)
chi_all(j,i,iparm)=chi(j,i)
! Restore sidechain parameters
do i=1,ntyp
do j=1,ntyp
- aa(j,i)=aa_all(j,i,iparm)
- bb(j,i)=bb_all(j,i,iparm)
+ aa_aq(j,i)=aa_all(j,i,iparm)
+ bb_aq(j,i)=bb_all(j,i,iparm)
r0(j,i)=r0_all(j,i,iparm)
sigma(j,i)=sigma_all(j,i,iparm)
chi(j,i)=chi_all(j,i,iparm)
+wturn3*eello_turn3 &
+wturn6*eello_turn6+wel_loc*eel_loc &
+edihcnstr+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
#else
etot=wsc*evdw+wscp*evdw2 &
+welec*(ees+evdw1) &
+wturn3*eello_turn3 &
+wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
#endif
#ifdef MPI
!--------------------------------------------------------------------------------
subroutine cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,*)
-#define DEBUG
+!#define DEBUG
#ifdef DEBUG
use geometry, only:int_from_cart1
use geometry_data, only:vbld,rad2deg,theta,phi,alph,omeg
call xdrffloat_(ixdrf, rtime, iret)
print *,"rtime",rtime," iret",iret !d
call xdrffloat_(ixdrf, rpotE, iret)
- write (iout,*) "rpotE",rpotE," iret",iret !d
+! write (iout,*) "rpotE",rpotE," iret",iret !d
call flush(iout)
call xdrffloat_(ixdrf, ruconst, iret)
call xdrffloat_(ixdrf, rt_bath, iret)
#else
call xdrffloat(ixdrf, rtime, iret)
call xdrffloat(ixdrf, rpotE, iret)
- write (iout,*) "rpotE",rpotE," iret",iret !d
+! write (iout,*) "rpotE",rpotE," iret",iret !d
call flush(iout)
call xdrffloat(ixdrf, ruconst, iret)
call xdrffloat(ixdrf, rt_bath, iret)
call xdrfint(ixdrf, jhpb(j), iret)
enddo
call xdrfint(ixdrf, nprop, iret)
- write (iout,*) "nprop",nprop !d
+! write (iout,*) "nprop",nprop !d
if (it.gt.0 .and. nprop.ne.nprop_prev) then
write (iout,*) "Warning previous nprop",nprop_prev,&
" current",nprop
#endif
if (iret.eq.0) exit
itraj=mod(it,totraj(iR,iparm))
-#define DEBUG
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "ii",ii," itraj",itraj," it",it
#endif
#ifdef DEBUG
write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,2*nres+2)
#endif
-#undef DEBUG
+!#undef DEBUG
if (iret.eq.0) exit
if (itmp .ne. nres + nct - nnt + 1) then
write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
endif
time=rtime
- write (iout,*) "calling slice" !d
+! write (iout,*) "calling slice" !d
call flush(iout) !d
islice=slice(nstep(itraj),time,is,ie,ts,te)
- write (iout,*) "islice",islice !d
+! write (iout,*) "islice",islice !d
call flush(iout) !d
do i=1,nres
" conformations stored so far, slice",islice
enddo
call flush(iout)
-#undef DEBUG
+!#undef DEBUG
return
end subroutine cxread
!--------------------------------------------------------------------------------
use names, only:ntyp1
use geometry_data
- use energy_data, only:itype,dsc
+ use energy_data, only:itype,dsc,molnum
use geometry, only:int_from_cart1
! use
! include "DIMENSIONS"
include "mpif.h"
integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
#endif
- integer :: j,k,l,ii,itj,iprint
+ integer :: j,k,l,ii,itj,iprint,mnum
if (.not. check_conf) then
conf_check=.true.
return
endif
call int_from_cart1(.false.)
do j=nnt+1,nct
- if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. &
+ mnum=molnum(j)
+ if (mnum.eq.5) cycle
+ if (itype(j-1,mnum).ne.ntyp1 .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
(vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
if (iprint.gt.0) &
write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
endif
enddo
do j=nnt,nct
- itj=itype(j)
- if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. &
+ mnum=molnum(j)
+ if (mnum.eq.5) cycle
+ itj=itype(j,mnum)
+ if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
(vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
if (iprint.gt.0) &
write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
! Read molecular data.
!
use energy_data
- use geometry_data, only:nres,deg2rad,c,dc
+ use geometry_data, only:nres,deg2rad,c,dc,nres_molec
use control_data, only:iscode
use control, only:rescode,setup_var,init_int_table
use geometry, only:alloc_geo_arrays
!el integer :: rescode
!el real(kind=8) :: x(maxvar)
character(len=320) :: controlcard !,ucase
- integer,dimension(nres) :: itype_pdb !(maxres)
- integer :: i,j,i1,i2,it1,it2
+ integer,dimension(nres,5) :: itype_pdb !(maxres)
+ integer :: i,j,i1,i2,it1,it2,mnum
real(kind=8) :: scalscp
!el logical :: seq_comp
+ write(iout,*) "START MOLREAD"
+ call flush(iout)
+ do i=1,5
+ nres_molec(i)=0
+ print *,"nres_molec, initial",nres_molec(i)
+ enddo
+
call card_concat(controlcard,.true.)
call reada(controlcard,'SCAL14',scal14,0.4d0)
call reada(controlcard,'SCALSCP',scalscp,1.0d0)
call reada(controlcard,'TEMP0',temp0,300.0d0) !el
call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
r0_corr=cutoff_corr-delt_corr
- call readi(controlcard,"NRES",nres,0)
- allocate(sequence(nres+1))
+ call readi(controlcard,"NRES",nres_molec(1),0)
+ call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
+ call readi(controlcard,"NRES_CAT",nres_molec(5),0)
+ nres=0
+ do i=1,5
+ nres=nres_molec(i)+nres
+ enddo
+! allocate(sequence(nres+1))
!el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
call alloc_geo_arrays
call alloc_ener_arrays
!----------------------------
allocate(c(3,2*nres+2))
allocate(dc(3,0:2*nres+2))
- allocate(itype(nres+2))
+ allocate(itype(nres+2,5))
allocate(itel(nres+2))
+ if (.not. allocated(molnum)) allocate(molnum(nres+2))
!
! Zero out tableis.
do i=1,2*nres+2
enddo
enddo
do i=1,nres+2
- itype(i)=0
+ molnum(i)=0
+ do j=1,5
+ itype(i,j)=0
+ enddo
itel(i)=0
enddo
!--------------------------
!
+ allocate(sequence(nres+1))
iscode=index(controlcard,"ONE_LETTER")
if (nres.le.0) then
write (iout,*) "Error: no residues in molecule"
write (iout,*) "Error: too many residues",nres,maxres
endif
write(iout,*) 'nres=',nres
+! HERE F**** CHANGE
! Read sequence of the protein
if (iscode.gt.0) then
read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
endif
! Convert sequence to numeric code
- do i=1,nres
- itype(i)=rescode(i,sequence(i),iscode)
+ do i=1,nres_molec(1)
+ mnum=1
+ molnum(i)=1
+ itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+ enddo
+ do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
+ mnum=2
+ molnum(i)=2
+ itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+ enddo
+ do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
+ mnum=5
+ molnum(i)=5
+ itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
enddo
+
write (iout,*) "Numeric code:"
- write (iout,'(20i4)') (itype(i),i=1,nres)
+ write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
do i=1,nres-1
+ mnum=molnum(i)
#ifdef PROCOR
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
+ if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
#else
- if (itype(i).eq.ntyp1) then
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
#endif
itel(i)=0
#ifdef PROCOR
- else if (iabs(itype(i+1)).ne.20) then
+ else if (iabs(itype(i+1,mnum)).ne.20) then
#else
- else if (iabs(itype(i)).ne.20) then
+ else if (iabs(itype(i,mnum)).ne.20) then
#endif
itel(i)=1
else
enddo
write (iout,*) "ITEL"
do i=1,nres-1
- write (iout,*) i,itype(i),itel(i)
+ mnum=molnum(i)
+ write (iout,*) i,itype(i,mnum),itel(i)
enddo
+ write(iout,*)
call read_bridge
if (with_dihed_constr) then
nnt=1
nct=nres
- if (itype(1).eq.ntyp1) nnt=2
- if (itype(nres).eq.ntyp1) nct=nct-1
+ if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
+ if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
write(iout,*) 'NNT=',NNT,' NCT=',NCT
call setup_var
call init_int_table
write (iout,'(20i4)') (iss(i),i=1,ns)
write (iout,'(/a/)') 'Pre-formed links are:'
do i=1,nss
+ mnum=molnum(i)
i1=ihpb(i)-nres
i2=jhpb(i)-nres
- it1=itype(i1)
- it2=itype(i2)
+ it1=itype(i1,mnum)
+ it2=itype(i2,mnum)
write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
- restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',&
+ restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
dhpb(i),ebr,forcon(i)
enddo
endif
character(len=16) :: key
integer :: iparm
!el real(kind=8) :: ip,mp
- real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,&
+ real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
res1
integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
- integer :: nlobi,iblock,maxinter,iscprol
+ integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
!
! Body
!
wvdwpp=ww(16)
wbond=ww(18)
wsccor=ww(19)
+ wcatcat=ww(44)
+ wcatprot=ww(43)
endif
!
weights(17)=wbond
weights(18)=0 !scal14 !
weights(21)=wsccor
+ weights(43)=wcatprot
+ weights(44)=wcatcat
+
! el--------
call card_concat(controlcard,.false.)
call reads(controlcard,"SCPPAR",scpname_t,scpname)
open (iscpp,file=scpname_t,status='old')
rewind(iscpp)
+ call getenv_loc('IONPAR',ionname)
+ open (iion,file=ionname,status='old')
+ rewind(iion)
write (iout,*) "Parameter set:",iparm
write (iout,*) "Energy-term weights:"
do i=1,n_eneW
endif
enddo
#else
- read (ibond,*) ijunk,vbldp0,akp,rjunk
+ read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
do i=1,ntyp
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
j=1,nbondterm(i))
'inertia','Pstok'
write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
do i=1,ntyp
- write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
+ write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
vbldsc0(1,i),aksc(1,i),abond0(1,i)
do j=2,nbondterm(i)
write (iout,'(13x,3f10.5)') &
enddo
enddo
endif
+ if (.not. allocated(msc)) allocate(msc(ntyp1,5))
+ if (.not. allocated(restok)) allocate(restok(ntyp1,5))
+
+ do i=1,ntyp_molec(5)
+ read(iion,*) msc(i,5),restok(i,5)
+ print *,msc(i,5),restok(i,5)
+ enddo
+ ip(5)=0.2
+! isc(5)=0.2
+ read (iion,*) ncatprotparm
+ allocate(catprm(ncatprotparm,4))
+ do k=1,4
+ read (iion,*) (catprm(i,k),i=1,ncatprotparm)
+ enddo
+ print *, catprm
+
+
!----------------------------------------------------
allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
+ allocate(epslip(ntyp,ntyp))
do i=1,ntyp
do j=1,ntyp
augm(i,j)=0.0D0
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,molnum(i)),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,molnum(i)),sigma0(i),rr0(i),&
i=1,ntyp)
endif
! goto 50
read (isidep,*)(sigii(i),i=1,ntyp)
read (isidep,*)(chip(i),i=1,ntyp)
read (isidep,*)(alp(i),i=1,ntyp)
+ do i=1,ntyp
+ read (isidep,*)(epslip(i,j),j=i,ntyp)
+ enddo
! For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
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,molnum(i)),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,molnum(i)),sigma0(i),rr0(i),&
sigii(i),chip(i),alp(i),i=1,ntyp)
endif
case default
! Calculate the "working" parameters of SC interactions.
!el from module energy - COMMON.INTERACT-------
- allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
do i=1,ntyp1
do j=1,ntyp1
- aa(i,j)=0.0D0
- bb(i,j)=0.0D0
+ aa_aq(i,j)=0.0D0
+ bb_aq(i,j)=0.0D0
+ aa_lip(i,j)=0.0d0
+ bb_lip(i,j)=0.0d0
chi(i,j)=0.0D0
sigma(i,j)=0.0D0
r0(i,j)=0.0D0
do i=2,ntyp
do j=1,i-1
eps(i,j)=eps(j,i)
+ epslip(i,j)=epslip(j,i)
enddo
enddo
do i=1,ntyp
epsij=eps(i,j)
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
- aa(i,j)=epsij*rrij*rrij
- bb(i,j)=-sigeps*epsij*rrij
- aa(j,i)=aa(i,j)
- bb(j,i)=bb(i,j)
+ aa_aq(i,j)=epsij*rrij*rrij
+ bb_aq(i,j)=-sigeps*epsij*rrij
+ aa_aq(j,i)=aa_aq(i,j)
+ bb_aq(j,i)=bb_aq(i,j)
+ epsijlip=epslip(i,j)
+ sigeps=dsign(1.0D0,epsijlip)
+ epsijlip=dabs(epsijlip)
+ aa_lip(i,j)=epsijlip*rrij*rrij
+ bb_lip(i,j)=-sigeps*epsijlip*rrij
+ aa_lip(j,i)=aa_lip(i,j)
+ bb_lip(j,i)=bb_lip(i,j)
if (ipot.gt.2) then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
endif
if (lprint) then
write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
- restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
+ restyp(i,molnum(i)),restyp(j,molnum(i)),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
!-----------------------------------------------------------------------------
subroutine read_general_data(*)
- use control_data, only:indpdb,symetr
+ use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions
use energy_data, only:distchainmax
+ use geometry_data, only:boxxsize,boxysize,boxzsize
! implicit none
! include "DIMENSIONS"
! include "DIMENSIONS.ZSCOPT"
call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
call readi(controlcard,"RESCALE",rescale_modeW,1)
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ ions=index(controlcard,"IONS").gt.0
+ call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+ call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
+ write(iout,*) "R_CUT_ELE=",r_cut_ele
check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
call readi(controlcard,'SYM',symetr,1)
use control_data, only:pdbref
use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
nstart_sup,nstart_seq,nperm,nres0
- use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
+ use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
use compare, only:seq_comp !,contact,elecont
use geometry, only:chainbuild,dist
use io_config, only:readpdb
character(len=4) :: sequence(nres)
!el integer rescode
!el real(kind=8) :: x(maxvar)
- integer :: itype_pdb(nres)
+ integer :: itype_pdb(nres,5)
!el logical seq_comp
- integer :: i,j,k,nres_pdb,iaux
+ integer :: i,j,k,nres_pdb,iaux,mnum
real(kind=8) :: ddsc !el,dist
integer :: kkk !,ilen
!el external ilen
return 1
34 continue
do i=1,nres
- itype_pdb(i)=itype(i)
+ mnum=molnum(i)
+ itype_pdb(i,mnum)=itype(i,mnum)
enddo
call readpdb
do i=1,nres
- iaux=itype_pdb(i)
- itype_pdb(i)=itype(i)
- itype(i)=iaux
+ iaux=itype_pdb(i,mnum)
+ itype_pdb(i,mnum)=itype(i,mnum)
+ itype(i,mnum)=iaux
enddo
close (ipdbin)
do kkk=1,nperm
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),&
+ if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
nsup)) then
do j=nnt+nsup-1,nnt,-1
do k=1,3
return 1
else
do i=0,nsup-(nct-nnt+1)
- if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
+ if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
nct-nnt+1)) &
then
nstart_sup=nstart_sup+i
enddo
enddo
do i=1,nres
+ mnum=molnum(i)
do j=1,3
dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
enddo
- if (itype(i).ne.10) then
+ if (itype(i,mnum).ne.10) then
ddsc = dist(i,nres+i)
do j=1,3
dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
use geometry_data, only:nres,c
- use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
+ use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
'D','E','F','G','H','I','J'/),shape(chainid))
integer,dimension(nres) :: ica !(maxres)
real(kind=8) :: temp,efree,etot,entropy,rmsdev
- integer :: ii,i,j,iti,ires,iatom,ichain
+ integer :: ii,i,j,iti,ires,iatom,ichain,mnum
write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
ii,temp,rmsdev
write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
ichain=1
ires=0
do i=nnt,nct
- iti=itype(i)
+ mnum=molnum(i)
+ iti=itype(i,mnum)
if (iti.eq.ntyp1) then
ichain=ichain+1
ires=0
ires=ires+1
iatom=iatom+1
ica(i)=iatom
- write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
+ write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,i),j=1,3)
if (iti.ne.10) then
iatom=iatom+1
- write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
+ write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,nres+i),j=1,3)
endif
endif
enddo
write (ipdb,'(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
+ mnum=molnum(i)
+ if (itype(i,mnum).eq.ntyp1) cycle
+ if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
write (ipdb,30) ica(i),ica(i+1)
- else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+ else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
write (ipdb,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,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
write (ipdb,30) ica(i),ica(i)+1
endif
enddo
- if (itype(nct).ne.10) then
+ if (itype(nct,molnum(nct)).ne.10) then
write (ipdb,30) ica(nct),ica(nct)+1
endif
do i=1,nss
esccor=enetb(21,i,iparm)
! edihcnstr=enetb(20,i,iparm)
edihcnstr=enetb(19,i,iparm)
+ ecationcation=enetb(42,i,iparm)
+ ecation_prot=enetb(41,i,iparm)
#ifdef DEBUG
write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),&
evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
+wturn3*eello_turn3 &
+wturn6*eello_turn6+wel_loc*eel_loc &
+edihcnstr+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
#else
etot=wsc*evdw+wscp*evdw2 &
+welec*(ees+evdw1) &
+wturn3*eello_turn3 &
+wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
#endif
#ifdef DEBUG
! edihcnstr=enetb(20,t,iparm)
edihcnstr=enetb(19,t,iparm)
edihcnstr=0.0d0
+ ecationcation=enetb(42,t,iparm)
+ ecation_prot=enetb(41,t,iparm)
+
do k=0,nGridT
betaT=startGridT+k*delta_T
temper=betaT
+ft(2)*wturn3*eello_turn3 &
+ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
+edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
- +wbond*estr
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
+ftprim(1)*wtor*etors+ &
ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
+ft(2)*wturn3*eello_turn3 &
+ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
- +wbond*estr
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
+ftprim(1)*wtor*etors+ &
ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
module wham_data
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
- integer,parameter :: max_eneW=21
+ integer,parameter :: max_eneW=49
integer,parameter :: maxQ=1
integer,parameter :: maxQ1=MaxQ+2
integer,parameter :: max_parm=1
projects / unres4.git / commitdiff