!
implicit none
!-----------------------------------------------------------------------------
+! Max. number of contacts per residue
+! integer :: maxconts
+!-----------------------------------------------------------------------------
+! Max. number of derivatives of virtual-bond and side-chain vectors in theta
+! or phi.
+! integer :: maxdim
+!-----------------------------------------------------------------------------
+! Max. number of SC contacts
+! integer :: maxcont
+!-----------------------------------------------------------------------------
! Max. number of variables
integer :: maxvar
!-----------------------------------------------------------------------------
-! Max number of torsional terms in SCCOR
- integer,parameter :: maxterm_sccor=6
+! Max number of torsional terms in SCCOR in control_data
+! integer,parameter :: maxterm_sccor=6
!-----------------------------------------------------------------------------
! Maximum number of SC local term fitting function coefficiants
integer,parameter :: maxsccoef=65
!-----------------------------------------------------------------------------
+! commom.calc common/calc/
+!-----------------------------------------------------------------------------
! commom.contacts
! common /contacts/
! Change 12/1/95 - common block CONTACTS1 included.
!
! Compute the side-chain and electrostatic interaction energy
!
- goto (101,102,103,104,105,106) ipot
+! goto (101,102,103,104,105,106) ipot
+ select case(ipot)
! Lennard-Jones potential.
- 101 call elj(evdw)
+! 101 call elj(evdw)
+ case (1)
+ call elj(evdw)
!d print '(a)','Exit ELJcall el'
- goto 107
+! goto 107
! Lennard-Jones-Kihara potential (shifted).
- 102 call eljk(evdw)
- goto 107
+! 102 call eljk(evdw)
+ case (2)
+ call eljk(evdw)
+! goto 107
! Berne-Pechukas potential (dilated LJ, angular dependence).
- 103 call ebp(evdw)
- goto 107
+! 103 call ebp(evdw)
+ case (3)
+ call ebp(evdw)
+! goto 107
! Gay-Berne potential (shifted LJ, angular dependence).
- 104 call egb(evdw)
- goto 107
+! 104 call egb(evdw)
+ case (4)
+ call egb(evdw)
+! goto 107
! Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
- 105 call egbv(evdw)
- goto 107
+! 105 call egbv(evdw)
+ case (5)
+ call egbv(evdw)
+! goto 107
! Soft-sphere potential
- 106 call e_softsphere(evdw)
+! 106 call e_softsphere(evdw)
+ case (6)
+ call e_softsphere(evdw)
!
! Calculate electrostatic (H-bonding) energy of the main chain.
!
- 107 continue
+! 107 continue
+ case default
+ write(iout,*)"Wrong ipot"
+! return
+! 50 continue
+ end select
+! continue
!mc
!mc Sep-06: egb takes care of dynamic ss bonds too
.or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
#endif
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+! write (iout,*) "ELEC calc"
else
ees=0.0d0
evdw1=0.0d0
! Calculate excluded-volume interaction energy between peptide groups
! and side chains.
!
+!elwrite(iout,*) "in etotal calc exc;luded",ipot
if (ipot.lt.6) then
if(wscp.gt.0d0) then
! write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
+!elwrite(iout,*) "in etotal before ebond",ipot
!
! Calculate the bond-stretching energy
!
call ebond(estr)
+!elwrite(iout,*) "in etotal afer ebond",ipot
+
!
! Calculate the disulfide-bridge and other energy and the contributions
! from other distance constraints.
- print *,'Calling EHPB'
+! print *,'Calling EHPB'
call edis(ehpb)
+!elwrite(iout,*) "in etotal afer edis",ipot
! print *,'EHPB exitted succesfully.'
!
! Calculate the virtual-bond-angle energy.
! Calculate the SC local energy.
!
call esc(escloc)
+!elwrite(iout,*) "in etotal afer esc",ipot
! print *,"Processor",myrank," computed USC"
!
! Calculate the virtual-bond torsional energy.
!
! 6/23/01 Calculate double-torsional energy
!
+!elwrite(iout,*) "in etotal",ipot
if (wtor_d.gt.0) then
call etor_d(etors_d)
else
ecorr6=0.0d0
eturn6=0.0d0
endif
+!elwrite(iout,*) "in etotal",ipot
if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
!d write (iout,*) "multibody_hb ecorr",ecorr
endif
+!elwrite(iout,*) "afeter multibody hb"
! print *,"Processor",myrank," computed Ucorr"
!
! If performing constraint dynamics, call the constraint energy
! after the equilibration time
if(usampl.and.totT.gt.eq_time) then
+!elwrite(iout,*) "afeter multibody hb"
call EconstrQ
+!elwrite(iout,*) "afeter multibody hb"
call Econstr_back
+!elwrite(iout,*) "afeter multibody hb"
else
Uconst=0.0d0
Uconst_back=0.0d0
endif
+!elwrite(iout,*) "after Econstr"
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#ifdef TIMING
time_sumene=time_sumene+MPI_Wtime()-time00
#endif
+!el call enerprint(energia)
+!elwrite(iout,*)"finish etotal"
return
end subroutine etotal
!-----------------------------------------------------------------------------
integer :: ierr
real(kind=8) :: time00
if (nfgtasks.gt.1 .and. reduce) then
-!el #define DEBUG
+
#ifdef DEBUG
write (iout,*) "energies before REDUCE"
call enerprint(energia)
#ifdef MPI
endif
#endif
-!el #undef DUBUG
+! call enerprint(energia)
call flush(iout)
return
end subroutine sum_energy
real(kind=8) :: kfac=2.4d0
real(kind=8) :: x,x2,x3,x4,x5,licznik=1.12692801104297249644
!el local variables
- real(kind=8) :: t_bath,facT,facT2,facT3,facT4,facT5
+ real(kind=8) :: t_bath,facT(6) !,facT2,facT3,facT4,facT5,facT6
+ real(kind=8) :: T0=3.0d2
integer :: ierror
! facT=temp0/t_bath
! facT=2*temp0/(t_bath+temp0)
if (rescale_mode.eq.0) then
- facT=1.0d0
- facT2=1.0d0
- facT3=1.0d0
- facT4=1.0d0
- facT5=1.0d0
+ facT(1)=1.0d0
+ facT(2)=1.0d0
+ facT(3)=1.0d0
+ facT(4)=1.0d0
+ facT(5)=1.0d0
+ facT(6)=1.0d0
else if (rescale_mode.eq.1) then
- facT=kfac/(kfac-1.0d0+t_bath/temp0)
- facT2=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2)
- facT3=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3)
- facT4=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4)
- facT5=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5)
+ facT(1)=kfac/(kfac-1.0d0+t_bath/temp0)
+ facT(2)=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2)
+ facT(3)=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3)
+ facT(4)=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4)
+ facT(5)=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5)
+#ifdef WHAM_RUN
+!#if defined(WHAM_RUN) || defined(CLUSTER)
+#if defined(FUNCTH)
+! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
+ facT(6)=(320.0+80.0*dtanh((t_bath-320.0)/80.0))/320.0
+#elif defined(FUNCT)
+ facT(6)=t_bath/T0
+#else
+ facT(6)=1.0d0
+#endif
+#endif
else if (rescale_mode.eq.2) then
x=t_bath/temp0
x2=x*x
x3=x2*x
x4=x3*x
x5=x4*x
- facT=licznik/dlog(dexp(x)+dexp(-x))
- facT2=licznik/dlog(dexp(x2)+dexp(-x2))
- facT3=licznik/dlog(dexp(x3)+dexp(-x3))
- facT4=licznik/dlog(dexp(x4)+dexp(-x4))
- facT5=licznik/dlog(dexp(x5)+dexp(-x5))
+ facT(1)=licznik/dlog(dexp(x)+dexp(-x))
+ facT(2)=licznik/dlog(dexp(x2)+dexp(-x2))
+ facT(3)=licznik/dlog(dexp(x3)+dexp(-x3))
+ facT(4)=licznik/dlog(dexp(x4)+dexp(-x4))
+ facT(5)=licznik/dlog(dexp(x5)+dexp(-x5))
+#ifdef WHAM_RUN
+!#if defined(WHAM_RUN) || defined(CLUSTER)
+#if defined(FUNCTH)
+ facT(6)=(320.0+80.0*dtanh((t_bath-320.0)/80.0))/320.0
+#elif defined(FUNCT)
+ facT(6)=t_bath/T0
+#else
+ facT(6)=1.0d0
+#endif
+#endif
else
write (iout,*) "Wrong RESCALE_MODE",rescale_mode
write (*,*) "Wrong RESCALE_MODE",rescale_mode
#endif
stop 555
endif
- welec=weights(3)*fact
- wcorr=weights(4)*fact3
- wcorr5=weights(5)*fact4
- wcorr6=weights(6)*fact5
- wel_loc=weights(7)*fact2
- wturn3=weights(8)*fact2
- wturn4=weights(9)*fact3
- wturn6=weights(10)*fact5
- wtor=weights(13)*fact
- wtor_d=weights(14)*fact2
- wsccor=weights(21)*fact
+ welec=weights(3)*fact(1)
+ wcorr=weights(4)*fact(3)
+ wcorr5=weights(5)*fact(4)
+ wcorr6=weights(6)*fact(5)
+ wel_loc=weights(7)*fact(2)
+ wturn3=weights(8)*fact(2)
+ wturn4=weights(9)*fact(3)
+ wturn6=weights(10)*fact(5)
+ wtor=weights(13)*fact(1)
+ wtor_d=weights(14)*fact(2)
+ wsccor=weights(21)*fact(1)
return
end subroutine rescale_weights
integer :: iint,itypi,itypi1,itypj
real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0ij
-
+ integer :: ii
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
evdw=evdw+evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
'evdw',i,j,evdwij,' ss'
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
ELSE
!el ind=ind+1
itypj=iabs(itype(j))
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,*) "j",j,dsc_inv(itypj),dscj_inv,&
+! 1.0d0/vbld(j+nres) !d
! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
! write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
-! write (iout,*) "j",j," dc_norm",
-! & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
+! write (iout,*) "j",j," dc_norm",& !d
+! dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
+! write(iout,*)"rrij ",rrij
+! write(iout,*)"xj yj zj ", xj, yj, zj
+! write(iout,*)"xi yi zi ", xi, yi, zi
+! write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
! Calculate angle-dependent terms of energy and contributions to their
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
rij_shift=1.0D0/rij-sig+sig0ij
+! write(iout,*)" rij_shift",rij_shift," rij",rij," sig",sig,&
+! "sig0ij",sig0ij
! for diagnostics; uncomment
! rij_shift=1.2*sig0ij
! I hate to put IF's in the loops, but here don't have another choice!!!!
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
-! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,&
-! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+! write(iout,*)"aa, bb ",aa(:,:),bb(:,:)
+! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d
+! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij
if (lprn) then
endif
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
- 'evdw',i,j,evdwij
+ 'evdw',i,j,evdwij !,"egb"
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
integer :: i,iti1,iti,k,l
real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2
-! allocate(Ug(2,2,nres)) !(2,2,maxres)
-! allocate(Ug2(2,2,nres)) !(2,2,maxres)
-! allocate(Ugder(2,2,nres)) !(2,2,maxres)
-! allocate(Ug2der(2,2,nres)) !(2,2,maxres)
-! allocate(obrot(2,nres)) !(2,maxres)
-! allocate(obrot2(2,nres)) !(2,maxres)
-! allocate(obrot_der(2,nres)) !(2,maxres)
-! allocate(obrot2_der(2,nres)) !(2,maxres)
-! allocate(costab2(nres)) !(maxres)
-! allocate(sintab2(nres)) !(maxres)
-! allocate(costab(nres)) !(maxres)
-! allocate(sintab(nres)) !(maxres)
-
-! allocate(Ub2(2,nres)) !(2,maxres)
-! allocate(Ctobr(2,nres)) !(2,maxres)
-! allocate(Dtobr2(2,nres)) !(2,maxres)
-! allocate(mu(2,nres)) !(2,maxres)
-! allocate(muder(2,nres)) !(2,maxres)
-! allocate(Ub2der(2,nres)) !(2,maxres)
-! allocate(Ctobrder(2,nres)) !(2,maxres)
-! allocate(Dtobr2der(2,nres)) !(2,maxres)
-
-! allocate(EUg(2,2,nres)) !(2,2,maxres)
-! allocate(CUg(2,2,nres)) !(2,2,maxres)
-! allocate(DUg(2,2,nres)) !(2,2,maxres)
-! allocate(DtUg2(2,2,nres)) !(2,2,maxres)
-! allocate(EUgder(2,2,nres)) !(2,2,maxres)
-! allocate(CUgder(2,2,nres)) !(2,2,maxres)
-! allocate(DUgder(2,2,nres)) !(2,2,maxres)
-! allocate(Dtug2der(2,2,nres)) !(2,2,maxres)
-
-! allocate(Ug2Db1t(2,nres)) !(2,maxres)
-! allocate(Ug2Db1tder(2,nres)) !(2,maxres)
-! allocate(CUgb2(2,nres)) !(2,maxres)
-! allocate(CUgb2der(2,nres)) !(2,maxres)
-
-! allocate(EUgC(2,2,nres)) !(2,2,maxres)
-! allocate(EUgCder(2,2,nres)) !(2,2,maxres)
-! allocate(EUgD(2,2,nres)) !(2,2,maxres)
-! allocate(EUgDder(2,2,nres)) !(2,2,maxres)
-! allocate(DtUg2EUg(2,2,nres)) !(2,2,maxres)
-! allocate(Ug2DtEUg(2,2,nres)) !(2,2,maxres)
-
-! allocate(Ug2DtEUgder(2,2,2,nres)) !(2,2,2,maxres)
-! allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres)
-
!
! Compute the virtual-bond-torsional-angle dependent quantities needed
! to calculate the el-loc multibody terms of various order.
!
+!AL el mu=0.0d0
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
enddo
-!d write (iout,*) 'mu ',mu(:,i-2)
+! if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2)
+! if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1)
+! if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2)
!d write (iout,*) 'mu1',mu1(:,i-2)
!d write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
enddo
endif
#if defined(MPI) && defined(PARMAT)
-!el #define DUBUG
#ifdef DEBUG
! if (fg_rank.eq.0) then
write (iout,*) "Arrays UG and UGDER before GATHER"
!d enddo
!d enddo
return
-!el #undef DUBUG
end subroutine set_matrices
!-----------------------------------------------------------------------------
subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
!d enddo
!d call check_vecgrad
!d stop
+! ees=0.0d0 !AS
+! evdw1=0.0d0
+! eel_loc=0.0d0
+! eello_turn3=0.0d0
+! eello_turn4=0.0d0
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ eel_loc=0.0d0
+ eello_turn3=0.0d0
+ eello_turn4=0.0d0
+!
+
if (icheckgrad.eq.1) then
+!el
+! do i=0,2*nres+2
+! dc_norm(1,i)=0.0d0
+! dc_norm(2,i)=0.0d0
+! dc_norm(3,i)=0.0d0
+! enddo
do i=1,nres-1
fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
do k=1,3
! Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+a33*muij(4)
-!d write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eelloc',i,j,eel_loc_ij
+! if (energy_dec) write (iout,*) "a22",a22," a23",a23," a32",a32," a33",a33
+! if (energy_dec) write (iout,*) "muij",muij
! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
eel_loc=eel_loc+eel_loc_ij
r0ij=2.20D0*rpp(iteli,itelj)
! r0ij=1.55D0*rpp(iteli,itelj)
call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
+!elwrite(iout,*) "num_conti",num_conti, "maxconts",maxconts
if (fcont.gt.0.0D0) then
num_conti=num_conti+1
if (num_conti.gt.maxconts) then
if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
#ifdef OSF
- phii=phi(i)
+ phii=phi(i)
if (phii.ne.phii) phii=150.0
#else
phii=phi(i)
endif
if (i.lt.nres .and. itype(i).ne.ntyp1) then
#ifdef OSF
- phii1=phi(i+1)
+ phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
phii1=pinorm(phii1)
z(1)=cos(phii1)
sumene1x,sumene2x,sumene3x,sumene4x,&
sumene1y,sumene2y,sumene3y,sumene4y,cossc,cossc1,&
cosfac2xx,sinfac2yy
-!el #define DEBUG
#ifdef DEBUG
real(kind=8) :: aincr,xxsave,sumenep,de_dxx_num,yysave,&
de_dyy_num,zzsave,de_dzz_num,costsave,sintsave,&
#ifdef DEBUG
write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
#endif
-!#undef DEBUG
!
!
cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
1 continue
enddo
-!el #undef DUBUG
return
end subroutine esc
!-----------------------------------------------------------------------------
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))
+ itori=itortyp(itype(i-2))
+ itori1=itortyp(itype(i-1))
phii=phi(i)
gloci=0.0D0
! Proline-Proline pair is a special case...
esccor_ii=0.0D0
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
+
! write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
phii=phi(i)
do intertyp=1,3 !intertyp
icall=0
call geom_to_var(nvar,x)
if (.not.split_ene) then
-write(iout,*) 'Calling CHECK_ECARTINT if'
call etotal(energia)
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
etot=energia(0)
!el call enerprint(energia)
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
call flush(iout)
write (iout,*) "enter cartgrad"
call flush(iout)
call cartgrad
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
write (iout,*) "exit cartgrad"
call flush(iout)
icall =1
do j=1,3
grad_s(j,0)=gcart(j,0)
enddo
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
do i=1,nres
do j=1,3
grad_s(j,i)=gcart(j,i)
enddo
enddo
else
-write(iout,*) 'Calling CHECK_ECARTIN else.'
!- split gradient check
call zerograd
call etotal_long(energia)
call zerograd
aincr=1.0D-7
print '(a)','Calling CHECK_INT.'
-write(iout,*) 'Calling CHECK_INT.'
nf=0
nfl=0
icg=1
call geom_to_var(nvar,x)
call var_to_geom(nvar,x)
call chainbuild
-write(iout,*) 'Calling CHECK_INT.'
icall=1
print *,'ICG=',ICG
call etotal(energia)
nfl=3
!d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
- write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
+!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
icall=1
do i=1,nvar
xi=x(i)
i,key,ii,gg(i),gana(i),&
100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
enddo
-write(iout,*) "jestesmy sobie w check eint!!"
return
end subroutine check_eint
!-----------------------------------------------------------------------------
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'evdw',i,j,evdwij
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,"egb_long"
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'evdw',i,j,evdwij
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,"egb_short"
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
! Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+a33*muij(4)
-!d write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eelloc',i,j,eel_loc_ij
-! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d
eel_loc=eel_loc+eel_loc_ij
! Partial derivatives in virtual-bond dihedral angles gamma
call sum_gradient
#ifdef TIMING
#endif
+!el write (iout,*) "After sum_gradient"
#ifdef DEBUG
+!el write (iout,*) "After sum_gradient"
do i=1,nres-1
write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
enddo
endif
+!elwrite (iout,*) "After sum_gradient"
#ifdef TIMING
time01=MPI_Wtime()
#endif
call intcartderiv
+!elwrite (iout,*) "After sum_gradient"
#ifdef TIMING
time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
#endif
endif
enddo
!alculate derivative of Tauangle
- do i=1,nres-1
- do j=1,3
- dc_norm2(j,i+nres)=-dc_norm(j,i+nres)
- enddo
- enddo
#ifdef PARINTDER
do i=itau_start,itau_end
#else
cost=dcos(theta(i))
cost1=dcos(omicron(2,i-1))
cosg=dcos(tauangle(1,i))
+!elwrite(iout,*) " vecpr5",i,nres
do j=1,3
+!elwrite(iout,*) " vecpreee",i,nres,j,i-2+nres
+!elwrite(iout,*) " vecpr5",dc_norm2(1,1)
dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
enddo
integer :: kstart,kend,lstart,lend,idummy
real(kind=8) :: delta=1.0d-7
!el local variables
- integer :: i,ii,j
+ integer :: i,ii,j
! real(kind=8) ::
! For the backbone
do i=0,nres-1
real(kind=8) :: deps,ssx0,ljx0
!-------END TESTING CODE
+ eij=0.0d0
i=resi
j=resj
endif
if (havebond) then
-#ifndef CLUST
-#ifndef WHAM
+!#ifndef CLUST
+!#ifndef WHAM
! if (dyn_ssbond_ij(i,j).eq.1.0d300) then
! write(iout,'(a15,f12.2,f8.1,2i5)')
! & "SSBOND_E_FORM",totT,t_bath,i,j
! endif
-#endif
-#endif
+!#endif
+!#endif
dyn_ssbond_ij(i,j)=eij
else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
dyn_ssbond_ij(i,j)=1.0d300
-#ifndef CLUST
-#ifndef WHAM
+!#ifndef CLUST
+!#ifndef WHAM
! write(iout,'(a15,f12.2,f8.1,2i5)')
! & "SSBOND_E_BREAK",totT,t_bath,i,j
-#endif
-#endif
+!#endif
+!#endif
endif
!-------TESTING CODE
- if (checkstop) then
+!el if (checkstop) then
if (jcheck.eq.0) write(iout,'(a,3f15.8,$)') &
"CHECKSTOP",rij,eij,ed
echeck(jcheck)=eij
- endif
+!el endif
enddo
if (checkstop) then
write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
! include 'COMMON.CHAIN'
! include 'COMMON.IOUNITS'
! include 'COMMON.SETUP'
-#ifndef CLUST
-#ifndef WHAM
! include 'COMMON.MD'
-#endif
-#endif
! Local variables
real(kind=8) :: emin
integer :: i,j,imin,ierr
!-----------------------------------------------------------------------------
#ifdef WHAM
subroutine read_ssHist
- implicit none
+! implicit none
! Includes
! include 'DIMENSIONS'
! include "DIMENSIONS.FREE"
maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond
!----------------------
! arrays in subroutine init_int_table
+!el#ifdef MPI
+!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)
+!el#endif
allocate(nint_gr(nres))
allocate(nscp_gr(nres))
allocate(ielstart(nres))
- allocate(ielend(nres)) !(maxres)
+ allocate(ielend(nres))
+!(maxres)
allocate(istart(nres,maxint_gr))
- allocate(iend(nres,maxint_gr)) !(maxres,maxint_gr)
+ allocate(iend(nres,maxint_gr))
+!(maxres,maxint_gr)
allocate(iscpstart(nres,maxint_gr))
- allocate(iscpend(nres,maxint_gr)) !(maxres,maxint_gr)
+ allocate(iscpend(nres,maxint_gr))
+!(maxres,maxint_gr)
allocate(ielstart_vdw(nres))
- allocate(ielend_vdw(nres)) !(maxres)
+ allocate(ielend_vdw(nres))
+!(maxres)
- allocate(lentyp(0:nfgtasks-1)) !(0:maxprocs-1)
+ allocate(lentyp(0:nfgtasks-1))
+!(0:maxprocs-1)
!----------------------
! commom.contacts
! common /contacts/
if(.not.allocated(icont_ref)) allocate(icont_ref(2,maxcont))
- allocate(icont(2,maxcont)) !(2,maxcont)
+ allocate(icont(2,maxcont))
+!(2,maxcont)
! common /contacts1/
- allocate(num_cont(0:nres+4)) !(maxres)
- allocate(jcont(maxconts,nres)) !(maxconts,maxres)
- allocate(facont(maxconts,nres)) !(maxconts,maxres)
- allocate(gacont(3,maxconts,nres)) !(3,maxconts,maxres)
+ allocate(num_cont(0:nres+4))
+!(maxres)
+ allocate(jcont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(facont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(gacont(3,maxconts,nres))
+!(3,maxconts,maxres)
! common /contacts_hb/
allocate(gacontp_hb1(3,maxconts,nres))
allocate(gacontp_hb2(3,maxconts,nres))
allocate(gacontm_hb2(3,maxconts,nres))
allocate(gacontm_hb3(3,maxconts,nres))
allocate(gacont_hbr(3,maxconts,nres))
- allocate(grij_hb_cont(3,maxconts,nres)) !(3,maxconts,maxres)
+ allocate(grij_hb_cont(3,maxconts,nres))
+!(3,maxconts,maxres)
allocate(facont_hb(maxconts,nres))
allocate(ees0p(maxconts,nres))
allocate(ees0m(maxconts,nres))
- allocate(d_cont(maxconts,nres)) !(maxconts,maxres)
- allocate(num_cont_hb(nres)) !(maxres)
- allocate(jcont_hb(maxconts,nres)) !(maxconts,maxres)
+ allocate(d_cont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(num_cont_hb(nres))
+!(maxres)
+ allocate(jcont_hb(maxconts,nres))
+!(maxconts,maxres)
! common /rotat/
allocate(Ug(2,2,nres))
allocate(Ugder(2,2,nres))
allocate(Ug2(2,2,nres))
- allocate(Ug2der(2,2,nres)) !(2,2,maxres)
+ allocate(Ug2der(2,2,nres))
+!(2,2,maxres)
allocate(obrot(2,nres))
allocate(obrot2(2,nres))
allocate(obrot_der(2,nres))
- allocate(obrot2_der(2,nres)) !(2,maxres)
+ allocate(obrot2_der(2,nres))
+!(2,maxres)
! common /precomp1/
allocate(mu(2,nres))
allocate(muder(2,nres))
allocate(Ub2(2,nres))
+ do i=1,nres
+ Ub2(1,i)=0.0d0
+ Ub2(2,i)=0.0d0
+ enddo
allocate(Ub2der(2,nres))
allocate(Ctobr(2,nres))
allocate(Ctobrder(2,nres))
allocate(Dtobr2(2,nres))
- allocate(Dtobr2der(2,nres)) !(2,maxres)
+ allocate(Dtobr2der(2,nres))
+!(2,maxres)
allocate(EUg(2,2,nres))
allocate(EUgder(2,2,nres))
allocate(CUg(2,2,nres))
allocate(DUg(2,2,nres))
allocate(Dugder(2,2,nres))
allocate(DtUg2(2,2,nres))
- allocate(DtUg2der(2,2,nres)) !(2,2,maxres)
+ allocate(DtUg2der(2,2,nres))
+!(2,2,maxres)
! common /precomp2/
allocate(Ug2Db1t(2,nres))
allocate(Ug2Db1tder(2,nres))
allocate(CUgb2(2,nres))
- allocate(CUgb2der(2,nres)) !(2,maxres)
+ allocate(CUgb2der(2,nres))
+!(2,maxres)
allocate(EUgC(2,2,nres))
allocate(EUgCder(2,2,nres))
allocate(EUgD(2,2,nres))
allocate(EUgDder(2,2,nres))
allocate(DtUg2EUg(2,2,nres))
- allocate(Ug2DtEUg(2,2,nres)) !(2,2,maxres)
+ allocate(Ug2DtEUg(2,2,nres))
+!(2,2,maxres)
allocate(Ug2DtEUgder(2,2,2,nres))
- allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres)
+ allocate(DtUg2EUgder(2,2,2,nres))
+!(2,2,2,maxres)
! common /rotat_old/
allocate(costab(nres))
allocate(sintab(nres))
allocate(costab2(nres))
- allocate(sintab2(nres)) !(maxres)
+ allocate(sintab2(nres))
+!(maxres)
! common /dipmat/
allocate(a_chuj(2,2,maxconts,nres))
!(2,2,maxconts,maxres)(maxconts=maxres/4)
allocate(ncont_sent(nres))
allocate(ncont_recv(nres))
- allocate(iat_sent(nres)) !(maxres)
+ allocate(iat_sent(nres))
+!(maxres)
allocate(iint_sent(4,nres,nres))
- allocate(iint_sent_local(4,nres,nres)) !(4,maxres,maxres)
+ allocate(iint_sent_local(4,nres,nres))
+!(4,maxres,maxres)
allocate(iturn3_sent(4,0:nres+4))
allocate(iturn4_sent(4,0:nres+4))
allocate(iturn3_sent_local(4,nres))
- allocate(iturn4_sent_local(4,nres)) !(4,maxres)
+ allocate(iturn4_sent_local(4,nres))
+!(4,maxres)
allocate(itask_cont_from(0:nfgtasks-1))
- allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
+ allocate(itask_cont_to(0:nfgtasks-1))
+!(0:max_fg_procs-1)
+
+
!----------------------
! commom.deriv;
! common /derivat/
allocate(dcdv(6,maxdim))
- allocate(dxdv(6,maxdim)) !(6,maxdim)
- allocate(dxds(6,nres)) !(6,maxres)
+ allocate(dxdv(6,maxdim))
+!(6,maxdim)
+ allocate(dxds(6,nres))
+!(6,maxres)
allocate(gradx(3,nres,0:2))
- allocate(gradc(3,nres,0:2)) !(3,maxres,2)
+ allocate(gradc(3,nres,0:2))
+!(3,maxres,2)
allocate(gvdwx(3,nres))
allocate(gvdwc(3,nres))
allocate(gelc(3,nres))
allocate(gcorr6_turn_long(3,nres))
allocate(gradxorr(3,nres))
allocate(gradcorr5(3,nres))
- allocate(gradcorr6(3,nres)) !(3,maxres)
+ allocate(gradcorr6(3,nres))
+!(3,maxres)
allocate(gloc(0:maxvar,0:2))
- allocate(gloc_x(0:maxvar,2)) !(maxvar,2)
+ allocate(gloc_x(0:maxvar,2))
+!(maxvar,2)
allocate(gel_loc(3,nres))
allocate(gel_loc_long(3,nres))
allocate(gcorr3_turn(3,nres))
allocate(gcorr4_turn(3,nres))
allocate(gcorr6_turn(3,nres))
allocate(gradb(3,nres))
- allocate(gradbx(3,nres)) !(3,maxres)
+ allocate(gradbx(3,nres))
+!(3,maxres)
allocate(gel_loc_loc(maxvar))
allocate(gel_loc_turn3(maxvar))
allocate(gel_loc_turn4(maxvar))
allocate(gel_loc_turn6(maxvar))
allocate(gcorr_loc(maxvar))
allocate(g_corr5_loc(maxvar))
- allocate(g_corr6_loc(maxvar)) !(maxvar)
+ allocate(g_corr6_loc(maxvar))
+!(maxvar)
allocate(gsccorc(3,nres))
- allocate(gsccorx(3,nres)) !(3,maxres)
- allocate(gsccor_loc(nres)) !(maxres)
- allocate(dtheta(3,2,nres)) !(3,2,maxres)
+ allocate(gsccorx(3,nres))
+!(3,maxres)
+ allocate(gsccor_loc(nres))
+!(maxres)
+ allocate(dtheta(3,2,nres))
+!(3,2,maxres)
allocate(gscloc(3,nres))
- allocate(gsclocx(3,nres)) !(3,maxres)
+ allocate(gsclocx(3,nres))
+!(3,maxres)
allocate(dphi(3,3,nres))
allocate(dalpha(3,3,nres))
- allocate(domega(3,3,nres)) !(3,3,maxres)
+ allocate(domega(3,3,nres))
+!(3,3,maxres)
! common /deriv_scloc/
allocate(dXX_C1tab(3,nres))
allocate(dYY_C1tab(3,nres))
allocate(dZZ_Ctab(3,nres))
allocate(dXX_XYZtab(3,nres))
allocate(dYY_XYZtab(3,nres))
- allocate(dZZ_XYZtab(3,nres)) !(3,maxres)
+ allocate(dZZ_XYZtab(3,nres))
+!(3,maxres)
! common /mpgrad/
allocate(jgrad_start(nres))
- allocate(jgrad_end(nres)) !(maxres)
+ allocate(jgrad_end(nres))
+!(maxres)
+!----------------------
! common /indices/
allocate(ibond_displ(0:nfgtasks-1))
allocate(iset_displ(0:nfgtasks-1))
allocate(iset_count(0:nfgtasks-1))
allocate(iint_count(0:nfgtasks-1))
- allocate(iint_displ(0:nfgtasks-1)) !(0:max_fg_procs-1)
+ allocate(iint_displ(0:nfgtasks-1))
+!(0:max_fg_procs-1)
!----------------------
! common.MD
! common /mdgrad/
allocate(gcart(3,0:nres))
- allocate(gxcart(3,0:nres)) !(3,0:MAXRES)
+ allocate(gxcart(3,0:nres))
+!(3,0:MAXRES)
allocate(gradcag(3,nres))
- allocate(gradxag(3,nres)) !(3,MAXRES)
+ allocate(gradxag(3,nres))
+!(3,MAXRES)
! common /back_constr/
!el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back)
allocate(dutheta(nres))
- allocate(dugamma(nres)) !(maxres)
+ allocate(dugamma(nres))
+!(maxres)
allocate(duscdiff(3,nres))
- allocate(duscdiffx(3,nres)) !(3,maxres)
+ allocate(duscdiffx(3,nres))
+!(3,maxres)
!el i io:read_fragments
! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20)
! allocate((:,:,:),allocatable :: ifrag_back !(3,maxfrag_back,maxprocs/20)
allocate(dUdconst(3,0:nres))
allocate(dUdxconst(3,0:nres))
allocate(dqwol(3,0:nres))
- allocate(dxqwol(3,0:nres)) !(3,0:MAXRES)
+ allocate(dxqwol(3,0:nres))
+!(3,0:MAXRES)
!----------------------
! common.sbridge
! common /sbridge/ in io_common: read_bridge
!el integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane
! common /dyn_ssbond/
! and side-chain vectors in theta or phi.
- allocate(dyn_ssbond_ij(0:nres+4,0:nres+4)) !(maxres,maxres)
+ allocate(dyn_ssbond_ij(0:nres+4,0:nres+4))
+!(maxres,maxres)
do i=1,nres
do j=i+1,nres
dyn_ssbond_ij(i,j)=1.0d300
enddo
if (nss.gt.0) then
- allocate(idssb(nss),jdssb(nss)) !(maxdim)
+ allocate(idssb(nss),jdssb(nss))
+!(maxdim)
endif
- allocate(dyn_ss_mask(nres)) !(maxres)
+ allocate(dyn_ss_mask(nres))
+!(maxres)
do i=1,nres
dyn_ss_mask(i)=.false.
enddo
! allocate(vlor2sccor(maxterm_sccor,20,20))
! allocate(vlor3sccor(maxterm_sccor,20,20)) !(maxterm_sccor,20,20)
!----------------
- allocate(gloc_sc(3,0:2*nres,0:10)) !(3,0:maxres2,10)maxres2=2*maxres
+ allocate(gloc_sc(3,0:2*nres,0:10))
+!(3,0:maxres2,10)maxres2=2*maxres
allocate(dcostau(3,3,3,2*nres))
allocate(dsintau(3,3,3,2*nres))
allocate(dtauangle(3,3,3,2*nres))
allocate(dcosomicron(3,3,3,2*nres))
- allocate(domicron(3,3,3,2*nres)) !(3,3,3,maxres2)maxres2=2*maxres
-!----------------------
-! common.scrot
-! Parameters of the SC rotamers (local) term
-! common/scrot/ in io_conf: parmread
-! allocate((:,:),allocatable :: sc_parmin !(maxsccoef,ntyp)
-!----------------------
-! common.torcnstr
-! common /torcnstr/
-!el in io_conf:molread
-! allocate((:),allocatable :: idih_constr,idih_nconstr !(maxdih_constr)
-! allocate((:),allocatable :: phi0,drange !(maxdih_constr)
-!----------------------
-! common.torsion
-! common/torsion/ in io_conf: parmread
-! allocate((:,:,:),allocatable :: v0 !(-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:,:),allocatable :: v1,v2 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:),allocatable :: vlor1 !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
-! allocate((:,:,:),allocatable :: vlor2,vlor3 !(maxlor,maxtor,maxtor)
-! allocate((:),allocatable :: itortyp !(-ntyp1:ntyp1)
-! allocate((:,:,:),allocatable :: nterm,nlor !(-maxtor:maxtor,-maxtor:maxtor,2)
-!
-! common /torsiond/ in io_conf: parmread
-! allocate((:,:,:,:,:,:),allocatable :: v1c,v1s
- !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:,:,:,:),allocatable :: v2c,v2s
- !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:,:),allocatable :: ntermd_1,ntermd_2
- !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
-! common/fourier/ in io_conf: parmread
-! allocate((:,:),allocatable :: b1,b2,&
-! b1tilde !(2,-maxtor:maxtor)
-! allocate((:,:,:),allocatable :: cc,dd,ee,&
-! ctilde,dtilde !(2,2,-maxtor:maxtor)
+ allocate(domicron(3,3,3,2*nres))
+!(3,3,3,maxres2)maxres2=2*maxres
!----------------------
! common.var
! common /restr/
- allocate(varall(maxvar)) !(maxvar)(maxvar=6*maxres)
+ allocate(varall(maxvar))
+!(maxvar)(maxvar=6*maxres)
allocate(mask_theta(nres))
allocate(mask_phi(nres))
- allocate(mask_side(nres)) !(maxres)
+ allocate(mask_side(nres))
+!(maxres)
!----------------------
! common.vectors
! common /vectors/
allocate(uy(3,nres))
- allocate(uz(3,nres)) !(3,maxres)
+ allocate(uz(3,nres))
+!(3,maxres)
allocate(uygrad(3,3,2,nres))
- allocate(uzgrad(3,3,2,nres)) !(3,3,2,maxres)
+ allocate(uzgrad(3,3,2,nres))
+!(3,3,2,maxres)
return
end subroutine alloc_ener_arrays
projects / unres4.git / blobdiff