ugamma_cache(i,ntwx_cache)=ugamma(i)
uscdiff_cache(i,ntwx_cache)=uscdiff(i)
enddo
-
+C print *,'przed returnbox'
+ call returnbox
do i=1,nres*2
do j=1,3
c_cache(j,i,ntwx_cache)=c(j,i)
if(me.eq.king) close(irest2)
return
end
-
+c------------------------------------------
+ subroutine returnbox
+ include 'DIMENSIONS'
+ include 'mpif.h'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ include 'COMMON.REMD'
+ include 'COMMON.SETUP'
+ include 'COMMON.MUCA'
+ include 'COMMON.HAIRPIN'
+ j=1
+ chain_beg=1
+C do i=1,nres
+C write(*,*) 'initial', i,j,c(j,i)
+C enddo
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
+ chain_end=i
+ if (allareout.eq.1) then
+ ireturnval=int(c(j,i)/boxxsize)
+ if (c(j,i).le.0) ireturnval=ireturnval-1
+ do k=chain_beg,chain_end
+ c(j,k)=c(j,k)-ireturnval*boxxsize
+ c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
+ enddo
+ endif
+ chain_beg=i+1
+ allareout=1
+ else
+ if (int(c(j,i)/boxxsize).eq.0) allareout=0
+ endif
+ enddo
+ if (allareout.eq.1) then
+ ireturnval=int(c(j,i)/boxxsize)
+ if (c(j,i).le.0) ireturnval=ireturnval-1
+ do k=chain_beg,nres
+ c(j,k)=c(j,k)-ireturnval*boxxsize
+ c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
+ enddo
+ endif
+C NO JUMP
+C do i=1,nres
+C write(*,*) 'befor no jump', i,j,c(j,i)
+C enddo
+ nojumpval=0
+ do i=2,nres
+ if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+ difference=abs(c(j,i-1)-c(j,i))
+C print *,'diff', difference
+ if (difference.gt.boxxsize/2.0) then
+ if (c(j,i-1).gt.c(j,i)) then
+ nojumpval=1
+ else
+ nojumpval=-1
+ endif
+ else
+ nojumpval=0
+ endif
+ endif
+ c(j,k)=c(j,k)+nojumpval*boxxsize
+ c(j,k+nres)=c(j,k+nres)+nojumpval*boxxsize
+ enddo
+ nojumpval=0
+ do i=2,nres
+ if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+ difference=abs(c(j,i-1)-c(j,i))
+ if (difference.gt.boxxsize/2.0) then
+ if (c(j,i-1).gt.c(j,i)) then
+ nojumpval=1
+ else
+ nojumpval=-1
+ endif
+ else
+ nojumpval=0
+ endif
+ endif
+ c(j,k)=c(j,k)+nojumpval*boxxsize
+ c(j,k+nres)=c(j,k+nres)+nojumpval*boxxsize
+ enddo
+
+ do i=1,nres
+ write(*,*) 'after no jump', i,j,c(j,i)
+ enddo
+
+C NOW Y dimension
+ j=2
+ chain_beg=1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
+ chain_end=i
+ if (allareout.eq.1) then
+ ireturnval=int(c(j,i)/boxysize)
+ if (c(j,i).le.0) ireturnval=ireturnval-1
+ do k=chain_beg,chain_end
+ c(j,k)=c(j,k)-ireturnval*boxysize
+ c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
+ enddo
+ endif
+ chain_beg=i+1
+ allareout=1
+ else
+ if (int(c(j,i)/boxysize).eq.0) allareout=0
+ endif
+ enddo
+ if (allareout.eq.1) then
+ ireturnval=int(c(j,i)/boxysize)
+ if (c(j,i).le.0) ireturnval=ireturnval-1
+ do k=chain_beg,nres
+ c(j,k)=c(j,k)-ireturnval*boxysize
+ c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
+ enddo
+ endif
+ nojumpval=0
+ do i=2,nres
+ if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+ difference=abs(c(j,i-1)-c(j,i))
+ if (difference.gt.boxysize/2.0) then
+ if (c(j,i-1).gt.c(j,i)) then
+ nojumpval=1
+ else
+ nojumpval=-1
+ endif
+ else
+ nojumpval=0
+ endif
+ endif
+ c(j,k)=c(j,k)+nojumpval*boxysize
+ c(j,k+nres)=c(j,k+nres)+nojumpval*boxysize
+ enddo
+ nojumpval=0
+ do i=2,nres
+ if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+ difference=abs(c(j,i-1)-c(j,i))
+ if (difference.gt.boxysize/2.0) then
+ if (c(j,i-1).gt.c(j,i)) then
+ nojumpval=1
+ else
+ nojumpval=-1
+ endif
+ else
+ nojumpval=0
+ endif
+ endif
+ c(j,k)=c(j,k)+nojumpval*boxysize
+ c(j,k+nres)=c(j,k+nres)+nojumpval*boxysize
+ enddo
+
+ j=3
+ chain_beg=1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
+ chain_end=i
+ if (allareout.eq.1) then
+ ireturnval=int(c(j,i)/boxysize)
+ if (c(j,i).le.0) ireturnval=ireturnval-1
+ do k=chain_beg,chain_end
+ c(j,k)=c(j,k)-ireturnval*boxzsize
+ c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
+ enddo
+ endif
+ chain_beg=i+1
+ allareout=1
+ else
+ if (int(c(j,i)/boxzsize).eq.0) allareout=0
+ endif
+ enddo
+ if (allareout.eq.1) then
+ ireturnval=int(c(j,i)/boxzsize)
+ if (c(j,i).le.0) ireturnval=ireturnval-1
+ do k=chain_beg,nres
+ c(j,k)=c(j,k)-ireturnval*boxzsize
+ c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
+ enddo
+ endif
+ nojumpval=0
+ do i=2,nres
+ if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+ difference=abs(c(j,i-1)-c(j,i))
+ if (difference.gt.(boxzsize/2.0)) then
+ if (c(j,i-1).gt.c(j,i)) then
+ nojumpval=1
+ else
+ nojumpval=-1
+ endif
+ else
+ nojumpval=0
+ endif
+ endif
+ c(j,k)=c(j,k)+nojumpval*boxzsize
+ c(j,k+nres)=c(j,k+nres)+nojumpval*boxzsize
+ enddo
+ nojumpval=0
+ do i=2,nres
+ if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+ difference=abs(c(j,i-1)-c(j,i))
+ if (difference.gt.boxzsize/2.0) then
+ if (c(j,i-1).gt.c(j,i)) then
+ nojumpval=1
+ else
+ nojumpval=-1
+ endif
+ else
+ nojumpval=0
+ endif
+ endif
+ c(j,k)=c(j,k)+nojumpval*boxzsize
+ c(j,k+nres)=c(j,k+nres)+nojumpval*boxzsize
+ enddo
+
+ return
+ end
C Max. number of AA residues
integer maxres
c parameter (maxres=250)
- parameter (maxres=100)
+ parameter (maxres=500)
C Appr. max. number of interaction sites
integer maxres2
parameter (maxres2=2*maxres)
C Max number of symetries
integer maxsym,maxperm
- parameter (maxsym=5,maxperm=120)
+ parameter (maxsym=12,maxperm=120)
C Max. number of variables
integer maxvar
parameter (maxvar=4*maxres)
integer MaxR,MaxT_h
integer MaxSlice
parameter (Max_Parm=1)
- parameter (MaxQ=1,MaxQ1=MaxQ+2)
- parameter(MaxR=1,MaxT_h=32)
+ parameter (MaxQ=4,MaxQ1=MaxQ+2)
+ parameter(MaxR=1,MaxT_h=36)
parameter(MaxSlice=40)
integer MaxN
parameter (MaxN=100)
C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 0 0 702
+C 0 0 703
subroutine cinfo
include 'COMMON.IOUNITS'
write(iout,*)'++++ Compile info ++++'
- write(iout,*)'Version 0.0 build 702'
- write(iout,*)'compiled Mon Dec 3 05:37:30 2012'
+ write(iout,*)'Version 0.0 build 703'
+ write(iout,*)'compiled Thu May 15 10:35:39 2014'
write(iout,*)'compiled by aks255@matrix.chem.cornell.edu'
write(iout,*)'OS name: Linux '
write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 '
enddo
enddo
enddo
+ write (iout,*) "indstart(me1),indend(me1)"
+ &,indstart(me1),indend(me1)
do i=indstart(me1),indend(me1)
#else
do iparm=1,nParmSet
enddo
do i=1,ntot
#endif
+
read(ientout,rec=i,err=101)
& ((csingle(l,k),l=1,3),k=1,nres),
& ((csingle(l,k+nres),l=1,3),k=nnt,nct),
& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
& wtor_d,wsccor,wbond
#endif
+C write (iout,*) "tuz przed energia"
call etotal(energia(0),fT)
+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)
endif
endif
endif
+C write (iout,*) "Czy tu dochodze"
potE(iii+1,iparm)=energia(0)
do k=1,21
enetb(k,iii+1,iparm)=energia(k)
goto 106
C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
105 call egbv(evdw,evdw_t)
+C write(iout,*) 'po elektostatyce'
C
C Calculate electrostatic (H-bonding) energy of the main chain.
C
- 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-C
+ 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+C write(iout,*) 'po eelec'
+
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
c
c Calculate the bond-stretching energy
c
+
call ebond(estr)
-c write (iout,*) "estr",estr
+C write (iout,*) "estr",estr
C
C Calculate the disulfide-bridge and other energy and the contributions
C from other distance constraints.
C Calculate the virtual-bond-angle energy.
C
call ebend(ebe)
-cd print *,'Bend energy finished.'
+C print *,'Bend energy finished.'
C
C Calculate the SC local energy.
C
call esc(escloc)
-cd print *,'SCLOC energy finished.'
+C print *,'SCLOC energy finished.'
C
C Calculate the virtual-bond torsional energy.
C
do i=1,nres
num_cont_hb(i)=0
enddo
-cd print '(a)','Enter EELEC'
-cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+C print '(a)','Enter EELEC'
+C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
do i=1,nres
gel_loc_loc(i)=0.0d0
gcorr_loc(i)=0.0d0
enddo
do i=iatel_s,iatel_e
+ if (i.eq.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+2).eq.ntyp1
& .or. itype(i-1).eq.ntyp1
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=0
-c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
+ if (j.eq.1) cycle
if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
& .or.itype(j+2).eq.ntyp1
& .or.itype(j-1).eq.ntyp1
c &'evdw1',i,j,evdwij
c &,iteli,itelj,aaa,evdw1
-c write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+ write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
c write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
c & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
c & 1.0D0/dsqrt(rrmij),evdwij,eesij,
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
endif
-
+ write (iout,'(a7,i5,4f7.3)')
+ & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
enddo
estr=0.5d0*AKP*estr+estr1
c
c write (iout,*) ithet_start,ithet_end
do i=ithet_start,ithet_end
C if (itype(i-1).eq.ntyp1) cycle
+ if (i.le.2) cycle
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
ichir21=isign(1,itype(i))
ichir22=isign(1,itype(i))
endif
+ if (i.eq.3) then
+ y(1)=0.0D0
+ y(2)=0.0D0
+ else
if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
y(1)=0.0D0
y(2)=0.0D0
endif
+ endif
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
& E_theta,E_tc)
endif
etheta=etheta+ethetai
+c write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
+c & 'ebend',i,ethetai,theta(i),itype(i)
c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i),
c & rad2deg*phii,rad2deg*phii1,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
etheta=0.0D0
c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
do i=ithet_start,ithet_end
+C if (i.eq.2) cycle
C if (itype(i-1).eq.ntyp1) cycle
+ if (i.le.2) cycle
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
if (iabs(itype(i+1)).eq.20) iblock=2
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
+ if (i.eq.3) then
+ phii=0.0d0
+ ityp1=nthetyp+1
+ do k=1,nsingle
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ else
if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
sinph1(k)=0.0d0
enddo
endif
+ endif
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
common /sccalc/ time11,time12,time112,theti,it,nlobit
delta=0.02d0*pi
escloc=0.0D0
-c write (iout,'(a)') 'ESC'
+C write (iout,*) 'ESC'
do i=loc_start,loc_end
it=itype(i)
if (it.eq.ntyp1) cycle
if (it.eq.10) goto 1
nlobit=nlob(iabs(it))
c print *,'i=',i,' it=',it,' nlobit=',nlobit
-c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
+C write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
theti=theta(i+1)-pipol
x(1)=dtan(theti)
x(2)=alph(i)
dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
enddo
dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
-c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c & esclocbi,ss,ssd
+ write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
+ & esclocbi,ss,ssd
escloci=ss*escloci+(1.0d0-ss)*esclocbi
c escloci=esclocbi
c write (iout,*) escloci
enddo
dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c & esclocbi,ss,ssd
+c & esclocbi,ss,ssd
escloci=ss*escloci+(1.0d0-ss)*esclocbi
-c write (iout,*) escloci
+C write (iout,*) 'i=',i, escloci
else
call enesc(x,escloci,dersc,ddummy,.false.)
endif
escloc=escloc+escloci
-c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+C write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+ write (iout,'(a6,i5,0pf7.3)')
+ & 'escloc',i,escloci
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& wscloc*dersc(1)
c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
+ if (i.le.2) cycle
if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
& .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
C if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
c lprn=.true.
etors_d=0.0D0
do i=iphi_start,iphi_end-1
+ if (i.le.3) cycle
C if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
C & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
& eps_scp(ntyp,2),rscp(ntyp,2),eps_orig(ntyp,ntyp)
c 12/5/03 modified 09/18/03 Bond stretching parameters.
double precision vbldp0,vbldsc0,akp,aksc,abond0,distchainmax
+ &,vbldpDUM
integer nbondterm
common /stretch/ vbldp0,vbldsc0(maxbondterm,ntyp),akp,
& aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp),
& distchainmax,nbondterm(ntyp)
+ &,vbldpDUM
+
#else
loc_start=2
loc_end=nres-1
- ithet_start=3
+ ithet_start=3
ithet_end=nres
iphi_start=nnt+3
iphi_end=nct