& tabpermchain,nchain ,npermchain,ireschain,iz_sc,nres_chomo
double precision c,cref,crefjlee,dc,xloc,xrot,dc_norm,t,r,prod,rt,
& rmssing,anatemp,chomo
- common /chain/ c(3,maxres2+2),dc(3,maxres2),xloc(3,maxres),
+ common /chain/ c(3,maxres2+2),dc(3,0:maxres2),xloc(3,maxres),
& xrot(3,maxres),dc_norm(3,maxres2),nres,nres0
common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres),
& rt(3,3,maxres)
double precision potE(MaxStr_Proc,Max_Parm),entfac(MaxStr_Proc),
- & q(MaxQ+2,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm)
+ & q(MaxQ+6,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm)
integer einicheck
common /energies/ potE,entfac,q,enetb,einicheck
integer MaxR,MaxT_h,maxHdim
integer MaxSlice
parameter (Max_Parm=5)
- parameter (MaxQ=4,MaxQ1=MaxQ+2)
+ parameter (MaxQ=4,MaxQ1=MaxQ+6)
parameter(MaxR=8,MaxT_h=36)
parameter(MaxSlice=40)
parameter(maxHdim=200)
implicit none
include "DIMENSIONS"
include "COMMON.IOUNITS"
+ include "COMMON.CONTROL"
integer nchain,nres,itype(nres),chain_border(2,maxchain),
& chain_length(nchain),itemp(maxchain),
& npermchain,tabpermchain(maxchain,maxperm),
nchain_group=nchain_group+1
iieq=1
iequiv(iieq,nchain_group)=i
+ if (symetr.eq.1) then
do j=i+1,nchain
if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
c k=0
iieq=iieq+1
iequiv(iieq,nchain_group)=j
enddo
+ endif
nequiv(nchain_group)=iieq
enddo
write(iout,*) "Number of equivalent chain groups:",nchain_group
c write (iout,*) "Before boxshift"
c call flush(iout)
c Box shift
- call oligomer
+c call oligomer
c write (iout,*) "After oligomer"
c call flush(iout)
do i=1,nres
& eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp,
& sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
- & appc(2,2),bppc(2,2)
+ & appc(2,2),bppc(2,2),epp_(2,2),rpp_(2,2)
double precision elcutoff,elecutoff_14
integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
double precision econt(maxcont)
+ double precision boxshift
*
* Load the constants of peptide bond - peptide bond interactions.
* Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
*
* as of 7/06/91.
*
-c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
-c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
+c data epp_ / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+ data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
data elpp6c /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
data elpp3c / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
& "Constants of electrostatic interaction energy expression."
do i=1,2
do j=1,2
- rri=rpp(i,j)**6
+ rri=rpp_(i,j)**6
appc(i,j)=epp(i,j)*rri*rri
bppc(i,j)=-2.0*epp(i,j)*rri
ael6c(i,j)=elpp6c(i,j)*4.2**6
xmedi=xi+0.5*dxi
ymedi=yi+0.5*dyi
zmedi=zi+0.5*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
+c write (iout,*) "i",xmedi,ymedi,zmedi
do 4 j=i+2,ien-1
jj=iperm(j,ipermmin)
ind=ind+1
xj=c(1,jj)+0.5*dxj
yj=c(2,jj)+0.5*dyj
zj=c(3,jj)+0.5*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+c write (iout,*) "j",xj,yj,zj
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
+c write (iout,*) "j",xj,yj,zj
rij=xj*xj+yj*yj+zj*zj
- sss=sscale(sqrt(rij))
- sssgrad=sscagrad(sqrt(rij))
rrmij=1.0/(xj*xj+yj*yj+zj*zj)
rmij=sqrt(rrmij)
r3ij=rrmij*rmij
endif
ees=ees+eesij
evdw=evdw+evdwij*sss
+c write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
4 continue
1 continue
if (lprint) then
double precision tt
integer snk_p(MaxR,MaxT_h,Max_parm)
logical lerr
+ integer ncont,icont(2,maxcont),isecstr(maxres)
character*256 bprotfile_temp
+ double precision totlength
call opentmp(islice,ientout,bprotfile_temp)
iii=0
ii=0
errmsg_count=0
c write (iout,*) "enecalc: nparmset ",nparmset
c write (iout,*) "enecalc: tormode ",tor_mode
+ write (iout,*) "ns",ns," dyn_ss",dyn_ss,(iss(i),i=1,ns)
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+ enddo
+ endif
+ write (iout,*) "dyn_ss_mask",(dyn_ss_mask(i),i=1,nres)
#ifdef MPI
do iparm=1,nParmSet
do ib=1,nT_h(iparm)
anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
q(nQ+1,iii+1)=rmsnat(iii+1,ipermin)
endif
- q(nQ+2,iii+1)=gyrate(iii+1)
+c write (iout,*) iii+1,q(nQ+3,iii+1),q(nQ+4,iii+1),q(nQ+5,iii+1)
c fT=T0*beta_h(ib,ipar)*1.987D-3
c ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
if (rescale_mode.eq.1) then
& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
& wtor_d,wsccor,wbond
#endif
-C write (iout,*) "tuz przed energia"
+c write (iout,*) "tuz przed energia"
call etotal(energia(0),fT)
-C write (iout,*) "tuz za energia"
+c write (iout,*) "tuz za energia"
#ifdef DEBUG
write (iout,*) "Conformation",i
write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
& ((c(l,k+nres),l=1,3),k=nnt,nct)
c call intout
- call pdbout(indstart(me1)+iii,
- & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
+c call pdbout(indstart(me1)+iii,
+c & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
call enerprint(energia(0),fT)
errmsg_count=errmsg_count+1
if (errmsg_count.gt.maxerrmsg_count)
endif
enddo ! iparm
+ q(nQ+2,iii+1)=gyrate(iii+1)
+c 8/28/2020 Adam - determine the fraction of secondary structures.
+ call elecont(.false.,ncont,icont,nnt,nct-1,1)
+ call secondary2(.false.,.false.,ncont,icont,isecstr)
+#ifdef DEBUG
+ write (iout,*) "secondary structure"
+ write (iout,'(80i1)') (isecstr(k),k=1,nres)
+#endif
+ q(nQ+3,iii+1)=0.0d0
+ q(nQ+4,iii+1)=0.0d0
+ q(nQ+5,iii+1)=0.0d0
+ totlength=0.0d0
+ do k=nnt,nct
+ if (itype(k).eq.ntyp1) cycle
+ totlength=totlength+1.0d0
+ l=isecstr(k)
+ q(nQ+3+l,iii+1)=q(nQ+3+l,iii+1)+1.0d0
+ enddo
+ q(nQ+3,iii+1)=q(nQ+3,iii+1)/totlength
+ q(nQ+4,iii+1)=q(nQ+4,iii+1)/totlength
+ q(nQ+5,iii+1)=q(nQ+5,iii+1)/totlength
+c write (iout,*) "iii",iii," nssbond",nssbond,nss
+c q(nQ+6,iii+1)=nssbond
+ q(nQ+6,iii+1)=nss
iii=iii+1
if (q(1,iii).le.0.0d0 .and. indpdb.gt.0)
call set_shield_fac2
endif
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-C write(iout,*) 'po eelec'
+c write(iout,*) 'po eelec eello_turn4',eello_turn4
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& call eturn4(i,eello_turn4)
+c write (iout,*) "i",i," eello_turn4",eello_turn4
#ifdef FOURBODY
num_cont_hb(i)=num_conti
#endif
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
& dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
& num_conti,j1,j2
+ double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
+ common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
j=i+3
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
include 'COMMON.HEADER'
include 'COMMON.SBRIDGE'
character*50 tytul
- character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/
+ character*1 chainid(32) /'A','B','C','D','E','F','G','H','I','J',
+ & 'K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
+ & '1','2','3','4','5','6'/
dimension ica(maxres)
write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')
& ii,temp,rmsdev
iatom=0
ichain=1
ires=0
+ iti_prev=0
do i=nnt,nct
iti=itype(i)
if (iti.eq.ntyp1) then
- ichain=ichain+1
ires=0
- write (ipdb,'(a)') 'TER'
+ if (iti_prev.ne.ntyp1) then
+ write (ipdb,'(a)') 'TER'
+ ichain=ichain+1
+ endif
else
ires=ires+1
iatom=iatom+1
ica(i)=iatom
write (ipdb,10) iatom,restyp(iti),chainid(ichain),
- & ires,(c(j,i),j=1,3)
+ & ires,(c(j,i),j=1,3),1.0d0
if (iti.ne.10) then
iatom=iatom+1
write (ipdb,20) iatom,restyp(iti),chainid(ichain),
- & ires,(c(j,nres+i),j=1,3)
+ & ires,(c(j,nres+i),j=1,3),1.0d0
endif
endif
+ iti_prev=iti
enddo
write (ipdb,'(a)') 'TER'
do i=nnt,nct-1
write (ipdb,30) ica(nct),ica(nct)+1
endif
do i=1,nss
- if (dyn_ss) then
- write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
- else
- write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
- endif
+ write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
enddo
write (ipdb,'(a6)') 'ENDMDL'
- 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
- 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+ 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
+ 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
30 FORMAT ('CONECT',8I5)
return
end
& link_end_peak
double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss
logical dyn_ss,dyn_ss_mask
+ integer nssbond
common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,Ht,
& dyn_ssbond_ij(max_cyst,max_cyst),
- & idssb(maxss),jdssb(maxss)
+ & idssb(maxss),jdssb(maxss),nssbond
common /dyn_ss_logic/
& dyn_ss,dyn_ss_mask(maxres)
include "COMMON.NAMES"
include "COMMON.TIME1"
include "COMMON.TORCNSTR"
+ include "COMMON.SETUP"
C
C The following is just to define auxiliary variables used in angle conversion
C
n_viol = 0
n_gviol = 0
n_map = 0
+ nfgtasks = 1
#ifndef SPLITELE
nprint_ene=nprint_ene-1
#endif
& ctemper(:ilen(ctemper))//"pdb"
endif
open(ipdb,file=pdbname)
- write (iout,*) "Before reading nlist",nlist
+c write (iout,*) "Before reading nlist",nlist
do i=1,nlist
read (ientout,rec=iperm(i))
& ((csingle(l,k),l=1,3),k=1,nres),
do i=1,ncont
i1=icont(1,i)
j1=icont(2,i)
- if (i1.ge.nstart_sup .and. i1.le.nend_sup
- & .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
+c if (i1.ge.nstart_sup .and. i1.le.nend_sup
+c & .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
cd write (iout,*) "parallel",i1,j1
if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
ii1=i1
& "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
endif
endif
- endif
+c endif
endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
enddo
double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
c-------END TESTING CODE
-
+ nssbond=0
i=resi
j=resj
ici=icys(i)
icj=icys(j)
if (ici.eq.0 .or. icj.eq.0) then
- write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
+ write (iout,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
& "Attempt to create",
& " a disulfide link between non-cysteine residues ",restyp(i),i,
& restyp(j),j
& ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
c-------END TESTING CODE
+c write (iout,'(2(a,i5),4(a,f7.2))') "resi",resi," resj",resj,
+c & " ljxm",ljxm," ljxs",ljxs," ssxm",ssxm," rij",rij
if (rij.gt.ljxm) then
havebond=.false.
ljd=rij-ljXs
& -2.0D0*alf12*eps3der+sigder*sigsq_om12
else if (rij.lt.ssxm) then
havebond=.true.
+ nssbond=nssbond+1
+c write (iout,*) "ssMD: nssbond",nssbond
ssd=rij-ssXs
eij=ssA*ssd*ssd+ssB*ssd+ssC
eij=eij*sss
eom2= 2*akth*deltat2+pom1-om1*pom2
eom12=pom2
else
+c nssbond=nssbond+1
omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
d_ssxm(1)=0.5D0*akct/ssA
cgrad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
enddo
return
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
-C include 'COMMON.SETUP'
#ifndef CLUST
#ifndef WHAM
C include 'COMMON.MD'
logical found
integer i_newnss(1024),displ(0:1024)
integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
+ nfgtasks=1
allnss=0
do i=1,ns-1
sumE_p(i,iparm)=0.0d0
sumEbis_p(i,iparm)=0.0d0
sumEsq_p(i,iparm)=0.0d0
- do j=1,nQ+2
+ do j=1,nQ+6
sumQ_p(j,i,iparm)=0.0d0
sumQsq_p(j,i,iparm)=0.0d0
sumEQ_p(j,i,iparm)=0.0d0
sumE(i,iparm)=0.0d0
sumEbis(i,iparm)=0.0d0
sumEsq(i,iparm)=0.0d0
- do j=1,nQ+2
+ do j=1,nQ+6
sumQ(j,i,iparm)=0.0d0
sumQsq(j,i,iparm)=0.0d0
sumEQ(j,i,iparm)=0.0d0
sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
- do j=1,nQ+2
+ do j=1,nQ+6
sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
sumE(k,iparm)=sumE(k,iparm)+etot*weight
sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
- do j=1,nQ+2
+ do j=1,nQ+6
sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
& sumW(i,iparm)
sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
& -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
- do j=1,nQ+2
+ do j=1,nQ+6
sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
& -sumQ(j,i,iparm)**2
& (startGridT+i*delta_T))+potEmin
write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
& sumW(i,iparm),sumE(i,iparm)
- write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+ write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+6)
write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
- & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+ & (sumQsq(j,i,iparm),j=1,nQ+6),(sumEQ(j,i,iparm),j=1,nQ+6)
write (iout,*)
write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
& sumW(i,iparm),sumE(i,iparm)
- write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+ write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+6)
write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
- & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+ & (sumQsq(j,i,iparm),j=1,nQ+6),(sumEQ(j,i,iparm),j=1,nQ+6)
write (34,*)
call flush(34)
enddo
projects / unres.git / commitdiff