c if (dyn_ss) call dyn_set_nss
c print *,"Processor",myrank," computed USCSC"
+c write (iout,*) "SCSC computed OK"
+c call flush_(iout)
#ifdef TIMING
time01=MPI_Wtime()
#endif
#ifdef TIMING
time_vec=time_vec+MPI_Wtime()-time01
#endif
+C Introduction of shielding effect first for each peptide group
+C the shielding factor is set this factor is describing how each
+C peptide group is shielded by side-chains
+C the matrix - shield_fac(i) the i index describe the ith between i and i+1
+C write (iout,*) "shield_mode",shield_mode
+ if (shield_mode.gt.0) then
+ call set_shield_fac
+ endif
c print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
#ifdef SPLITELE
c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
c & eello_turn4)
endif
+c write (iout,*) "eelec computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed UELEC"
C
C Calculate excluded-volume interaction energy between peptide groups
c write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
+c write (iout,*) "escp computed OK"
+c call flush_(iout)
c
c Calculate the bond-stretching energy
c
call ebond(estr)
+c write (iout,*) "ebond computed OK"
+c call flush_(iout)
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
if (wang.gt.0d0) then
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
else
ebe=0
+ ethetacnstr=0
endif
+c write (iout,*) "ebend computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed UB"
C
C Calculate the SC local energy.
C
C print *,"TU DOCHODZE?"
call esc(escloc)
+c write (iout,*) "esc computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed USC"
C
C Calculate the virtual-bond torsional energy.
etors=0
edihcnstr=0
endif
+c write (iout,*) "etor computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed Utor"
C
C 6/23/01 Calculate double-torsional energy
else
etors_d=0
endif
+c write (iout,*) "etor_d computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed Utord"
C
C 21/5/07 Calculate local sicdechain correlation energy
else
esccor=0.0d0
endif
+c write (iout,*) "eback_sc_corr computed OK"
+c call flush_(iout)
C print *,"PRZED MULIt"
c print *,"Processor",myrank," computed Usccorr"
C
n_corr1=0
if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0
& .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then
+c write (iout,*) "Calling multibody_eello"
+c call flush_(iout)
call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
-cd write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
-cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+c write(iout,*)
+c & 'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
+c & " ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+c call flush_(iout)
else
ecorr=0.0d0
ecorr5=0.0d0
eturn6=0.0d0
endif
if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
+c write (iout,*) "Calling multibody_gb_ecorr"
+c call flush_(iout)
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
-cd write (iout,*) "multibody_hb ecorr",ecorr
+c write (iout,*) "Exited multibody_hb ecorr",ecorr
+c call flush_(iout)
endif
+c write (iout,*) "multibody computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed Ucorr"
C
C If performing constraint dynamics, call the constraint energy
if (wliptran.gt.0) then
call Eliptransfer(eliptran)
endif
-C print *,"za lipidami"
+c write (iout,*) "lipid energy computed OK"
+c call flush_(iout)
if (AFMlog.gt.0) then
call AFMforce(Eafmforce)
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
endif
+c write (iout,*) "AFMforce computed OK"
+c call flush_(iout)
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(21)=esccor
energia(22)=eliptran
energia(23)=Eafmforce
+ energia(24)=ethetacnstr
c Here are the energies showed per procesor if the are more processors
c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
+c write (iout,*) "sum energy OK"
+c call flush_(iout)
if (dyn_ss) call dyn_set_nss
+c write (iout,*) "Exiting energy"
+c call flush_(iout)
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
time_sumene=time_sumene+MPI_Wtime()-time00
esccor=energia(21)
eliptran=energia(22)
Eafmforce=energia(23)
+ ethetacnstr=energia(24)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
+ & +ethetacnstr
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
& +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
& +Eafmforce
+ & +ethetacnstr
#endif
energia(0)=etot
c detecting NaNQ
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
enddo
enddo
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
enddo
enddo
do i=-1,nct
do j=1,3
#ifdef SPLITELE
+C print *,gradbufc(1,13)
+C print *,welec*gelc(1,13)
+C print *,wel_loc*gel_loc(1,13)
+C print *,0.5d0*(wscp*gvdwc_scpp(1,13))
+C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13)
+C print *,wel_loc*gel_loc_long(1,13)
+C print *,gradafm(1,13),"AFM"
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
& +wscloc*gscloc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
+ & +welec*gshieldc_loc(j,i)
+
+
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& +wscloc*gscloc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
+ & +welec*gshieldc_loc(j,i)
+
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
& +wliptran*gliptranx(j,i)
+ & +welec*gshieldx(j,i)
enddo
enddo
#ifdef DEBUG
esccor=energia(21)
eliptran=energia(22)
Eafmforce=energia(23)
+ ethetacnstr=energia(24)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
- & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
- & edihcnstr,ebr*nss,
- & Uconst,eliptran,wliptran,Eafmforce,etot
+ & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
+ & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
+ & etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot
+ & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
+ & etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
& 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
include 'COMMON.SBRIDGE'
logical lprn
integer xshift,yshift,zshift
+
evdw=0.0D0
ccccc energy_dec=.false.
C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+
+c write(iout,*) "PRZED ZWYKLE", evdwij
call dyn_ssbond_ene(i,j,evdwij)
+c write(iout,*) "PO ZWYKLE", evdwij
+
evdw=evdw+evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
& 'evdw',i,j,evdwij,' ss'
+C triple bond artifac removal
+ do k=j+1,iend(i,iint)
+C search over all next residues
+ if (dyn_ss_mask(k)) then
+C check if they are cysteins
+C write(iout,*) 'k=',k
+
+c write(iout,*) "PRZED TRI", evdwij
+ evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
+c if(evdwij_przed_tri.ne.evdwij) then
+c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+c endif
+
+c write(iout,*) "PO TRI", evdwij
+C call the energy function that removes the artifical triple disulfide
+C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
ELSE
ind=ind+1
itypj=iabs(itype(j))
c write (iout,*) 'theta=', theta(i-1)
enddo
#else
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itortyp(itype(i-2))
+ else
+ iti=ntortyp+1
+ endif
+c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itortyp(itype(i-1))
+ else
+ iti1=ntortyp+1
+ endif
b1(1,i-2)=b(3,iti)
b1(2,i-2)=b(5,iti)
b2(1,i-2)=b(2,iti)
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
+C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
c write (iout,*) 'mu ',mu(:,i-2),i-2
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
if (i.le.1) cycle
C write(iout,*) "tu jest i",i
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+4).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes by Ana
& .or. itype(i+2).eq.ntyp1
- & .or. itype(i+3).eq.ntyp1
- & .or. itype(i-1).eq.ntyp1
- & .or. itype(i+4).eq.ntyp1
- & ) cycle
+ & .or. itype(i+3).eq.ntyp1) cycle
+ if(i.gt.1)then
+ if(itype(i-1).eq.ntyp1)cycle
+ end if
+ if(i.LT.nres-3)then
+ if (itype(i+4).eq.ntyp1) cycle
+ end if
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
do i=iturn4_start,iturn4_end
if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+5).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes suggested by Ana
& .or. itype(i+3).eq.ntyp1
& .or. itype(i+4).eq.ntyp1
& .or. itype(i+5).eq.ntyp1
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
+CTU KURWA
do i=iatel_s,iatel_e
+C do i=75,75
if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+2).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes by Ana
& .or. itype(i+2).eq.ntyp1
& .or. itype(i-1).eq.ntyp1
& ) cycle
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
+C I TU KURWA
do j=ielstart(i),ielend(i)
+C do j=16,17
C write (iout,*) i,j
if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((j+2).gt.nres)
+ & .or.((j-1).le.0)
+C end of changes by Ana
& .or.itype(j+2).eq.ntyp1
& .or.itype(j-1).eq.ntyp1
&) cycle
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
include 'COMMON.SPLITELE'
+ include 'COMMON.SHIELD'
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
el2=fac4*fac
C MARYSIA
- eesij=(el1+el2)
+C eesij=(el1+el2)
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
+ if (shield_mode.gt.0) then
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ el1=el1*fac_shield(i)*fac_shield(j)
+ el2=el2*fac_shield(i)*fac_shield(j)
+ eesij=(el1+el2)
ees=ees+eesij
+ else
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ eesij=(el1+el2)
+ ees=ees+eesij
+ endif
evdw1=evdw1+evdwij*sss
cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
+
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
ggg(1)=facel*xj
ggg(2)=facel*yj
ggg(3)=facel*zj
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C if (iresshield.gt.i) then
+C do ishi=i+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C else
+C do ishi=iresshield,i
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C endif
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C if (iresshield.gt.j) then
+C do ishi=j+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C
+C enddo
+C else
+C do ishi=iresshield,j
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C enddo
+C endif
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,i)+
+ & grad_shield(k,i)*eesij/fac_shield(i)
+ gshieldc(k,j)=gshieldc(k,j)+
+ & grad_shield(k,j)*eesij/fac_shield(j)
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & grad_shield(k,i)*eesij/fac_shield(i)
+ gshieldc(k,j-1)=gshieldc(k,j-1)+
+ & grad_shield(k,j)*eesij/fac_shield(j)
+
+ enddo
+ endif
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
c gelc(k,j)=gelc(k,j)+ghalf
c enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
+C print *,"before", gelc_long(1,i), gelc_long(1,j)
do k=1,3
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+C & +grad_shield(k,j)*eesij/fac_shield(j)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+C & +grad_shield(k,i)*eesij/fac_shield(i)
+C gelc_long(k,i-1)=gelc_long(k,i-1)
+C & +grad_shield(k,i)*eesij/fac_shield(i)
+C gelc_long(k,j-1)=gelc_long(k,j-1)
+C & +grad_shield(k,j)*eesij/fac_shield(j)
enddo
+C print *,"bafter", gelc_long(1,i), gelc_long(1,j)
+
*
* Loop over residues i+1 thru j-1.
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
ggg(1)=fac*xj
+C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j)
ggg(2)=fac*yj
+C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j)
ggg(3)=fac*zj
+C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j)
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
cd & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
+ & fac_shield(i)*fac_shield(j)
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
+C print *,"before22", gelc_long(1,i), gelc_long(1,j)
do k=1,3
gelc(k,i)=gelc(k,i)
- & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+ & *fac_shield(i)*fac_shield(j)
gelc(k,j)=gelc(k,j)
- & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+ & *fac_shield(i)*fac_shield(j)
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+C print *,"before33", gelc_long(1,i), gelc_long(1,j)
+
C MARYSIA
c endif !sscale
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
& +a33*muij(4)
c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
c & ' eel_loc_ij',eel_loc_ij
-c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
+C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
C Calculate patrial derivative for theta angle
#ifdef NEWCORR
geel_loc_ij=a22*gmuij1(1)
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
dimension ggg(3)
ehpb=0.0D0
+ do i=1,3
+ ggg(i)=0.0d0
+ enddo
+C write (iout,*) ,"link_end",link_end,constr_dist
cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
cd write(iout,*)'link_start=',link_start,' link_end=',link_end
if (link_end.eq.0) return
C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
if (.not.dyn_ss .and. i.le.nss) then
C 15/02/13 CC dynamic SSbond - additional check
- if (ii.gt.nres
- & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+ if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+ & iabs(itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
endif
cd write (iout,*) "eij",eij
+cd & ' waga=',waga,' fac=',fac
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (constr_dist.eq.11) then
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+ & ehpb,fordepth(i),dd
+ else
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
else
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
+ if (constr_dist.eq.11) then
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+ & ehpb,fordepth(i),dd
+ else
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
waga=forcon(i)
C Calculate the contribution to energy.
ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
fac=waga*rdis/dd
-cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-cd & ' waga=',waga,' fac=',fac
+ endif
+ endif
do j=1,3
ggg(j)=fac*(c(j,jj)-c(j,ii))
enddo
enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
return
end
C--------------------------------------------------------------------------
end
#ifdef CRYST_THETA
C--------------------------------------------------------------------------
- subroutine ebend(etheta)
+ subroutine ebend(etheta,ethetacnstr)
C
C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
C angles gamma and its derivatives in consecutive thetas and gammas.
include 'COMMON.NAMES'
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
+ include 'COMMON.TORCNSTR'
common /calcthet/ term1,term2,termm,diffak,ratak,
& ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
& delthe0,sig0inv,sigtc,sigsqtc,delthec,it
etheta=0.0D0
c write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
+c write (iout,*) "ebend: i=",i
+c call flush_(iout)
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.
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
+c write (iout,*) "Exit loop"
+c call flush_(iout)
+ ethetacnstr=0.0d0
+c write (iout,*) ithetaconstr_start,ithetaconstr_end,"TU"
+c call flush_(iout)
+ do i=max0(ithetaconstr_start,1),ithetaconstr_end
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+ & i,itheta,rad2deg*thetiii,
+ & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+ & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+ & gloc(itheta+nphi-2,icg)
+ endif
+ enddo
+c write (iout,*) "Exit ebend"
+c call flush_(iout)
+
C Ufff.... We've done all this!!!
return
end
end
#else
C--------------------------------------------------------------------------
- subroutine ebend(etheta)
+ subroutine ebend(etheta,ethetacnstr)
C
C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
C angles gamma and its derivatives in consecutive thetas and gammas.
include 'COMMON.NAMES'
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
+ include 'COMMON.TORCNSTR'
double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
& cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
& sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+ if (i.le.2) cycle
c print *,i,itype(i-1),itype(i),itype(i-2)
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
-C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
-
+C print *,i,theta(i)
if (iabs(itype(i+1)).eq.20) iblock=2
if (iabs(itype(i+1)).ne.20) iblock=1
dethetai=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
+C print *,ethetai
+ 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)
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
do k=1,nsingle
+ ityp1=ithetyp((itype(i-2)))
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
+ endif
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp((itype(i)))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
enddo
write(iout,*) "ethetai",ethetai
endif
+C print *,ethetai
do m=1,ntheterm2
do k=1,nsingle
aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
& ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
& ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
& eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
+C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
enddo
enddo
+C print *,"cosph1", (cosph1(k), k=1,nsingle)
+C print *,"cosph2", (cosph2(k), k=1,nsingle)
+C print *,"sinph1", (sinph1(k), k=1,nsingle)
+C print *,"sinph2", (sinph2(k), k=1,nsingle)
if (lprn)
& write(iout,*) "ethetai",ethetai
+C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
do m=1,ntheterm3
do k=2,ndouble
do l=1,k-1
enddo
10 continue
c lprn1=.true.
+C print *,ethetai
if (lprn1)
& write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
& i,theta(i)*rad2deg,phii*rad2deg,
etheta=etheta+ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
- gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
+ enddo
+C now constrains
+ ethetacnstr=0.0d0
+C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=max0(ithetaconstr_start,1),ithetaconstr_end
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+ & i,itheta,rad2deg*thetiii,
+ & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+ & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+ & gloc(itheta+nphi-2,icg)
+ endif
enddo
+
return
end
#endif
difi=phii-phi0(i)
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
endif
! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
-c lprn=.true.
+c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
C ANY TWO ARE DUMMY ATOMS in row CYCLE
difi=pinorm(phii-phi0(i))
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else
difi=0.0
endif
-cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
-cd & rad2deg*phi0(i), rad2deg*drange(i),
-cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc",
+ & i,itori,rad2deg*phii,
+ & rad2deg*phi0(i), rad2deg*drange(i),
+ & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg)
+ endif
enddo
cd write (iout,*) 'edihcnstr',edihcnstr
return
if (lprn) then
write (iout,'(a)') 'Contact function values before RECEIVE:'
do i=nnt,nct-2
- write (iout,'(2i3,50(1x,i2,f5.2))')
+ write (iout,'(2i3,50(1x,i3,f5.2))')
& i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
& j=1,num_cont_hb(i))
enddo
jp1=iabs(j1)
c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
c & ' jj=',jj,' kk=',kk
+c call flush(iout)
if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0
& .or. j.lt.0 .and. j1.gt.0) .and.
& (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then
enddo ! kk
do kk=1,num_conti
j1=jcont_hb(kk,i)
-c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
-c & ' jj=',jj,' kk=',kk
+c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+c & ' jj=',jj,' kk=',kk
+c call flush(iout)
if (j1.eq.j+1) then
C Contacts I-J and (I+1)-J occur simultaneously.
C The system loses extra energy.
include 'COMMON.CONTACTS'
include 'COMMON.CHAIN'
include 'COMMON.CONTROL'
+ include 'COMMON.TORSION'
double precision gx(3),gx1(3)
integer num_cont_hb_old(maxres)
logical lprn,ldone
C Set lprn=.true. for debugging
lprn=.false.
eturn6=0.0d0
+c write (iout,*) "MULTIBODY_EELLO"
+c call flush(iout)
#ifdef MPI
do i=1,nres
num_cont_hb_old(i)=num_cont_hb(i)
if (lprn) then
write (iout,'(a)') 'Contact function values before RECEIVE:'
do i=nnt,nct-2
- write (iout,'(2i3,50(1x,i2,f5.2))')
+ write (iout,'(2i3,50(1x,i3,f5.2))')
& i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
& j=1,num_cont_hb(i))
enddo
+ call flush(iout)
endif
- call flush(iout)
do i=1,ntask_cont_from
ncont_recv(i)=0
enddo
if (lprn) then
write (iout,'(a)') 'Contact function values:'
do i=nnt,nct-2
- write (iout,'(2i3,50(1x,i2,5f6.3))')
+ write (iout,'(2i3,50(1x,i3,5f6.3))')
& i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i),
& ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i))
enddo
+ write (iout,*) "itortyp"
+ do i=1,nres
+ write (iout,*) i,itype(i),itortyp(itype(i))
+ enddo
+ call flush(iout)
endif
ecorr=0.0D0
ecorr5=0.0d0
do kk=1,num_conti1
j1=jcont_hb(kk,i1)
jp1=iabs(j1)
-c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
-c & ' jj=',jj,' kk=',kk
+ if (lprn) then
+ write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+ & ' jj=',jj,' kk=',kk
+ call flush(iout)
+ endif
c if (j1.eq.j+1 .or. j1.eq.j-1) then
if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0
& .or. j.lt.0 .and. j1.gt.0) .and.
c write (iout,'(i5,3f10.5)')
c & iii,(gradcorr5(jjj,iii),jjj=1,3)
c enddo
+c write (iout,*) "ecorr4"
+c call flush(iout)
+c write (iout,*) "eello5:",i,jp,i+1,jp1,jj,kk,
+c & itype(jp),itype(i+1),itype(jp1),
+c & itortyp(itype(jp)),itortyp(itype(i+1)),itortyp(itype(jp1))
+c call flush(iout)
if (wcorr5.gt.0.0d0)
& ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk)
c write (iout,*) "gradcorr5 after eello5"
2 'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk)
cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6
cd write(2,*)'ijkl',i,jp,i+1,jp1
+c write (iout,*) "ecorr5"
+c call flush(iout)
if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3
& .or. wturn6.eq.0.0d0))then
cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk)
if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
1 'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk)
+c write (iout,*) "ecorr6"
+c call flush(iout)
cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
cd & 'ecorr6=',ecorr6
cd write (iout,'(4e15.5)') sred_geom,
if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
1 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk)
cd write (2,*) 'multibody_eello:eturn6',eturn6
+c write (iout,*) "ecorr4"
+c call flush(iout)
endif
ENDIF
1111 continue
endif
+ if (energy_dec) call flush(iout)
enddo ! kk
enddo ! jj
enddo ! i
cd write (iout,*)
cd & 'EELLO5: Contacts have occurred for peptide groups',i,j,
cd & ' and',k,l
- itk=itortyp(itype(k))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+c itk=itortyp(itype(k))
+c itl=itortyp(itype(l))
+c itj=itortyp(itype(j))
eello5_1=0.0d0
eello5_2=0.0d0
eello5_3=0.0d0
c goto 1112
c1111 continue
C Contribution from graph II
- call transpose2(EE(1,1,itk),auxmat(1,1))
+ call transpose2(EE(1,1,k),auxmat(1,1))
call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
cd goto 1112
C Contribution from graph IV
cd1110 continue
- call transpose2(EE(1,1,itl),auxmat(1,1))
+ call transpose2(EE(1,1,l),auxmat(1,1))
call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
cd goto 1112
C Contribution from graph IV
1110 continue
- call transpose2(EE(1,1,itj),auxmat(1,1))
+ call transpose2(EE(1,1,j),auxmat(1,1))
call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k)
cgrad ghalf=0.5d0*ggg2(ll)
cd ghalf=0.0d0
- gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2)
+ gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2)
gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2)
- gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2)
+ gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2)
gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2)
gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl
gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl
C i i C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- itk=itortyp(itype(k))
+c itk=itortyp(itype(k))
s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
C
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
C energy moment and not to the cluster cumulant.
- iti=itortyp(itype(i))
- if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
- else
- itj1=ntortyp
- endif
- itk=itortyp(itype(k))
- itk1=itortyp(itype(k+1))
- if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
- else
- itl1=ntortyp
- endif
+c iti=itortyp(itype(i))
+c if (j.lt.nres-1) then
+c itj1=itortyp(itype(j+1))
+c else
+c itj1=ntortyp
+c endif
+c itk=itortyp(itype(k))
+c itk1=itortyp(itype(k+1))
+c if (l.lt.nres-1) then
+c itl1=itortyp(itype(l+1))
+c else
+c itl1=ntortyp
+c endif
#ifdef MOMENT
s1=dip(4,jj,i)*dip(4,kk,k)
#endif
s2=0.5d0*scalar2(b1(1,k),auxvec(1))
call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
- call transpose2(EE(1,1,itk),auxmat(1,1))
+ call transpose2(EE(1,1,k),auxmat(1,1))
call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
vv(1)=pizda(1,1)+pizda(2,2)
vv(2)=pizda(2,1)-pizda(1,2)
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
C energy moment and not to the cluster cumulant.
cd write (2,*) 'eello_graph4: wturn6',wturn6
- iti=itortyp(itype(i))
- itj=itortyp(itype(j))
- if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
- else
- itj1=ntortyp
- endif
- itk=itortyp(itype(k))
- if (k.lt.nres-1) then
- itk1=itortyp(itype(k+1))
- else
- itk1=ntortyp
- endif
- itl=itortyp(itype(l))
- if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
- else
- itl1=ntortyp
- endif
+c iti=itortyp(itype(i))
+c itj=itortyp(itype(j))
+c if (j.lt.nres-1) then
+c itj1=itortyp(itype(j+1))
+c else
+c itj1=ntortyp
+c endif
+c itk=itortyp(itype(k))
+c if (k.lt.nres-1) then
+c itk1=itortyp(itype(k+1))
+c else
+c itk1=ntortyp
+c endif
+c itl=itortyp(itype(l))
+c if (l.lt.nres-1) then
+c itl1=itortyp(itype(l+1))
+c else
+c itl1=ntortyp
+c endif
cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l
cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,
cd & ' itl',itl,' itl1',itl1
C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
return
end
+C-----------------------------------------------------------
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+ double precision div77_81/0.974996043d0/,
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=dsqrt(dist_pep_side)
+ dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
+C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist
+ & *(2.0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5
+C remember for the final gradient multiply sh_frac_dist_grad(j)
+C for side_chain by factor -2 !
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C if ((i.eq.3).and.(k.eq.2)) then
+C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C & ,"TU"
+C endif
+
+C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k))
+ long=long_r_sidechain(itype(k))
+ costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
+C now costhet_grad
+C costhet=0.0d0
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+C remember for the final gradient multiply costhet_grad(j)
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+
+ cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+ enddo
+
+ VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
+ & /VSolvSphere_div
+C now the gradient...
+C grad_shield is gradient of Calfa for peptide groups
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)
+C gradient po costhet
+ &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+ &-scale_fac_dist*(cosphi_grad_long(j))
+ &/(1.0-cosphi) )*div77_81
+ &*VofOverlap
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*(-2.0d0)
+ & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+ & +scale_fac_dist*(cosphi_grad_long(j))
+ & *2.0d0/(1.0-cosphi))
+ & *div77_81*VofOverlap
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*cosphi_grad_loc(j)
+ & *2.0d0/(1.0-cosphi)
+ & *div77_81*VofOverlap
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*div77_81+div4_81
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end