real(kind=8),dimension(:,:,:),allocatable :: gacont !(3,maxconts,maxres)
integer,dimension(:),allocatable :: ishield_list
integer,dimension(:,:),allocatable :: shield_list
+ real(kind=8),dimension(:),allocatable :: enetube,enecavtube
!
! 12/26/95 - H-bonding contacts
! common /contacts_hb/
gshieldc_ec,gshieldc_loc_ec,gshieldx_t3, &
gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, &
gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,&
- grad_shield,gg_tube,gg_tube_sc !(3,maxres)
+ grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres)
! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
integer :: n_corr,n_corr1,ierror
real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
- real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube
+ real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
+ Eafmforce,ethetacnstr
real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6
#ifdef MPI
#endif
!
! Compute the side-chain and electrostatic interaction energy
- print *, "Before EVDW"
+! print *, "Before EVDW"
! goto (101,102,103,104,105,106) ipot
select case(ipot)
! Lennard-Jones potential.
! print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
#ifdef SPLITELE
- print *,"after ipot if", ipot
+! print *,"after ipot if", ipot
if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or. &
wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 &
.or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
! Calculate the virtual-bond-angle energy.
!
if (wang.gt.0d0) then
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
else
ebe=0
endif
else
eliptran=0.0d0
endif
+ if (fg_rank.eq.0) then
+ if (AFMlog.gt.0) then
+ call AFMforce(Eafmforce)
+ else if (selfguide.gt.0) then
+ call AFMvel(Eafmforce)
+ endif
+ endif
if (tubemode.eq.1) then
call calctube(etube)
else if (tubemode.eq.2) then
energia(20)=Uconst+Uconst_back
energia(21)=esccor
energia(22)=eliptran
+ energia(23)=Eafmforce
+ energia(24)=ethetacnstr
energia(25)=etube
! Here are the energies showed per procesor if the are more processors
! per molecule then we sum it up in sum_energy subroutine
real(kind=8) :: evdw,evdw2,evdw2_14,ees,evdw1,ecorr,ecorr5,ecorr6
real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc
real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot, &
- eliptran,etube
+ eliptran,etube, Eafmforce,ethetacnstr
integer :: i
#ifdef MPI
integer :: ierr
Uconst=energia(20)
esccor=energia(21)
eliptran=energia(22)
+ Eafmforce=energia(23)
+ ethetacnstr=energia(24)
etube=energia(25)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 &
+wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
+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+wtube*etube
+ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+ +Eafmforce+ethetacnstr
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
+wang*ebe+wtor*etors+wscloc*escloc &
+wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
+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+wtube*etube
+ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+ +Eafmforce+ethetacnstr
+
#endif
energia(0)=etot
! detecting NaNQ
real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ecorr,ecorr5,ecorr6,eel_loc
real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc
real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,&
- etube
+ etube,ethetacnstr,Eafmforce
etot=energia(0)
evdw=energia(1)
Uconst=energia(20)
esccor=energia(21)
eliptran=energia(22)
+ Eafmforce=energia(23)
+ ethetacnstr=energia(24)
etube=energia(25)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
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,etube,wtube,etot
+ edihcnstr,ethetacnstr,ebr*nss,&
+ Uconst,eliptran,wliptran,Eafmforce,etube,wtube,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)'/&
+ 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ &
'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
'ETOT= ',1pE16.6,' (total)')
#else
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,etube,wtube,etot
+ ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, &
+ etube,wtube,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)'/ &
+ 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ &
'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
'ETOT= ',1pE16.6,' (total)')
#endif
sslipi=0.0d0
ssgradlipi=0.0
endif
- print *, sslipi,ssgradlipi
+! print *, sslipi,ssgradlipi
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
'evdw',i,j,evdwij,' ss'
! if (energy_dec) write (iout,*) &
! 'evdw',i,j,evdwij,' ss'
+ 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
!el ind=ind+1
itypj=iabs(itype(j))
! write (iout,*) 'i',i,' fac',fac
enddo
endif
- print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4, &
- wturn6
+! print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4, &
+! wturn6
if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
.or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or. &
wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
ehpb=ehpb+2*eij
!d write (iout,*) "eij",eij
endif
+ 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
-! Calculate the distance between the two points and its difference from the
-! target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
-! Get the force constant corresponding to this distance.
- waga=forcon(i)
-! Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
-!
-! Evaluate gradient.
-!
- fac=waga*rdis/dd
-!d print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-!d & ' waga=',waga,' fac=',fac
- do j=1,3
- ggg(j)=fac*(c(j,jj)-c(j,ii))
- enddo
-!d print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
-! If this is a SC-SC distance, we need to calculate the contributions to the
-! Cartesian gradient in the SC vectors (ghpbx).
- if (iii.lt.ii) then
+ 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
+ endif
+ endif
+
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+!cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
+!C If this is a SC-SC distance, we need to calculate the contributions to the
+!C Cartesian gradient in the SC vectors (ghpbx).
+ if (iii.lt.ii) then
do j=1,3
ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
- endif
-!grad do j=iii,jjj-1
-!grad do k=1,3
-!grad ghpbc(k,j)=ghpbc(k,j)+ggg(k)
-!grad enddo
-!grad enddo
- do k=1,3
- ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
- ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
- enddo
+ endif
+!cgrad do j=iii,jjj-1
+!cgrad do k=1,3
+!cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+!cgrad enddo
+!cgrad enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+
return
end subroutine edis
!-----------------------------------------------------------------------------
usumsqder=usumsqder+ud(j)*uprod2
enddo
estr=estr+uprod/usum
+ if (energy_dec) write (iout,*) &
+ "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
+ AKSC(1,iti),AKSC(1,iti)*diff*diff
do j=1,3
gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
enddo
end subroutine theteng
#else
!-----------------------------------------------------------------------------
- subroutine ebend(etheta)
+ subroutine ebend(etheta,ethetacnstr)
!
! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
! angles gamma and its derivatives in consecutive thetas and gammas.
!el local variables
integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
- real(kind=8) :: aux,etheta,ccl,ssl,scl,csl
+ real(kind=8) :: aux,etheta,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+ real(kind=8) :: difi,thetiii
+ integer itheta
etheta=0.0D0
do i=ithet_start,ithet_end
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=wang*dethetai
enddo
+!-----------thete constrains
+! if (tor_mode.ne.2) then
+ ethetacnstr=0.0d0
+!C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=ithetaconstr_start,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
+! endif
+
return
end subroutine ebend
#endif
wturn6*gcorr6_turn_long(j,i)+ &
wstrain*ghpbc(j,i) &
+wliptran*gliptranc(j,i) &
+ +gradafm(j,i) &
+welec*gshieldc(j,i) &
+wcorr*gshieldc_ec(j,i) &
+wturn3*gshieldc_t3(j,i)&
wturn6*gcorr6_turn_long(j,i)+ &
wstrain*ghpbc(j,i) &
+wliptran*gliptranc(j,i) &
+ +gradafm(j,i) &
+welec*gshieldc(j,i)&
+wcorr*gshieldc_ec(j,i) &
+wturn4*gshieldc_t4(j,i) &
wsccor*gsccorc(j,i) &
+wscloc*gscloc(j,i) &
+wliptran*gliptranc(j,i) &
+ +gradafm(j,i) &
+welec*gshieldc(j,i) &
+welec*gshieldc_loc(j,i) &
+wcorr*gshieldc_ec(j,i) &
wturn6*gcorr6_turn(j,i)+ &
wsccor*gsccorc(j,i) &
+wscloc*gscloc(j,i) &
+ +gradafm(j,i) &
+wliptran*gliptranc(j,i) &
+welec*gshieldc(j,i) &
+welec*gshieldc_loc(j,) &
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
- call dyn_ssbond_ene(i,j,evdwij)
- evdw=evdw+evdwij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
- 'evdw',i,j,evdwij,' ss'
+! call dyn_ssbond_ene(i,j,evdwij)
+! 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'
+! 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
!el ind=ind+1
itypj=itype(j)
evdw=evdw+evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
'evdw',i,j,evdwij,' ss'
+ 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
+
! if (energy_dec) write (iout,*) &
! 'evdw',i,j,evdwij,' ss'
ELSE
!el local variables
integer :: i,nres6
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,esccor,etors_d,etors
- real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr
+ real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr
nres6=6*nres
! write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
!
! Calculate the virtual-bond-angle energy.
!
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
!
! Calculate the SC local energy.
!
endif
return
end function gnmr1prim
-!-----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+ real(kind=8) function rlornmr1(y,ymin,ymax,sigma)
+ real(kind=8) y,ymin,ymax,sigma
+ real(kind=8) wykl /4.0d0/
+ if (y.lt.ymin) then
+ rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl)
+ else if (y.gt.ymax) then
+ rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl)
+ else
+ rlornmr1=0.0d0
+ endif
+ return
+ end function rlornmr1
+!------------------------------------------------------------------------------
+ real(kind=8) function rlornmr1prim(y,ymin,ymax,sigma)
+ real(kind=8) y,ymin,ymax,sigma
+ real(kind=8) wykl /4.0d0/
+ if (y.lt.ymin) then
+ rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ &
+ ((ymin-y)**wykl+sigma**wykl)**2
+ else if (y.gt.ymax) then
+ rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ &
+ ((y-ymax)**wykl+sigma**wykl)**2
+ else
+ rlornmr1prim=0.0d0
+ endif
+ return
+ end function rlornmr1prim
+
real(kind=8) function harmonic(y,ymax)
! implicit none
real(kind=8) :: y,ymax
gshieldc_loc_ll(j,i)=0.0d0
gg_tube(j,i)=0.0d0
gg_tube_sc(j,i)=0.0d0
+ gradafm(j,i)=0.0d0
do intertyp=1,3
gloc_sc(intertyp,i,icg)=0.0d0
enddo
return
end subroutine dyn_ssbond_ene
+!--------------------------------------------------------------------------
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
+! implicit none
+! Includes
+ use calc_data
+ use comm_sschecks
+! include 'DIMENSIONS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ use MD_data
+! include 'COMMON.MD'
+! use MD, only: totT,t_bath
+#endif
+#endif
+ double precision h_base
+ external h_base
+
+!c Input arguments
+ integer resi,resj,resk,m,itypi,itypj,itypk
+
+!c Output arguments
+ double precision eij,eij1,eij2,eij3
+
+!c Local variables
+ logical havebond
+!c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+ double precision xik,yik,zik,xjk,yjk,zjk,dxk,dyk,dzk
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ eij=0.0
+ if (dtriss.eq.0) return
+ i=resi
+ j=resj
+ k=resk
+!C write(iout,*) resi,resj,resk
+ itypi=itype(i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ itypj=itype(j)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+ itypk=itype(k)
+ xk=c(1,nres+k)
+ yk=c(2,nres+k)
+ zk=c(3,nres+k)
+
+ dxk=dc_norm(1,nres+k)
+ dyk=dc_norm(2,nres+k)
+ dzk=dc_norm(3,nres+k)
+ dscj_inv=vbld_inv(k+nres)
+ xij=xj-xi
+ xik=xk-xi
+ xjk=xk-xj
+ yij=yj-yi
+ yik=yk-yi
+ yjk=yk-yj
+ zij=zj-zi
+ zik=zk-zi
+ zjk=zk-zj
+ rrij=(xij*xij+yij*yij+zij*zij)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+ rrik=(xik*xik+yik*yik+zik*zik)
+ rik=dsqrt(rrik)
+ rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+ rjk=dsqrt(rrjk)
+!C there are three combination of distances for each trisulfide bonds
+!C The first case the ith atom is the center
+!C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+!C distance y is second distance the a,b,c,d are parameters derived for
+!C this problem d parameter was set as a penalty currenlty set to 1.
+ if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
+ eij1=0.0d0
+ else
+ eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
+ endif
+!C second case jth atom is center
+ if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
+ eij2=0.0d0
+ else
+ eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
+ endif
+!C the third case kth atom is the center
+ if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
+ eij3=0.0d0
+ else
+ eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
+ endif
+!C eij2=0.0
+!C eij3=0.0
+!C eij1=0.0
+ eij=eij1+eij2+eij3
+!C write(iout,*)i,j,k,eij
+!C The energy penalty calculated now time for the gradient part
+!C derivative over rij
+ fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
+ gg(1)=xij*fac/rij
+ gg(2)=yij*fac/rij
+ gg(3)=zij*fac/rij
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,j)=gvdwx(m,j)+gg(m)
+ enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+!C now derivative over rik
+ fac=-eij1**2/dtriss* &
+ (-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xik*fac/rik
+ gg(2)=yik*fac/rik
+ gg(3)=zik*fac/rik
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+!C now derivative over rjk
+ fac=-eij2**2/dtriss* &
+ (-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- &
+ eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xjk*fac/rjk
+ gg(2)=yjk*fac/rjk
+ gg(3)=zjk*fac/rjk
+ do m=1,3
+ gvdwx(m,j)=gvdwx(m,j)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+ return
+ end subroutine triple_ssbond_ene
+
+
+
!-----------------------------------------------------------------------------
real(kind=8) function h_base(x,deriv)
! A smooth function going 0->1 in range [0,1]
diff=newnss-nss
!mc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
-
+! print *,newnss,nss,maxdim
do i=1,nss
found=.false.
+! print *,newnss
do j=1,newnss
+!! print *,j
if (idssb(i).eq.newihpb(j) .and. &
jdssb(i).eq.newjhpb(j)) found=.true.
enddo
#ifndef CLUST
#ifndef WHAM
+! write(iout,*) "found",found,i,j
if (.not.found.and.fg_rank.eq.0) &
write(iout,'(a15,f12.2,f8.1,2i5)') &
"SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
do i=1,newnss
found=.false.
do j=1,nss
+! print *,i,j
if (newihpb(i).eq.idssb(j) .and. &
newjhpb(i).eq.jdssb(j)) found=.true.
enddo
#ifndef CLUST
#ifndef WHAM
+! write(iout,*) "found",found,i,j
if (.not.found.and.fg_rank.eq.0) &
write(iout,'(a15,f12.2,f8.1,2i5)') &
"SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
integer :: i
eliptran=0.0
- print *, "I am in eliptran"
+! print *, "I am in eliptran"
do i=ilip_start,ilip_end
!C do i=1,1
if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))&
!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
!C simple Kihara potential
subroutine calctube(Etube)
- real(kind=8) :: vectube(3),enetube(nres*2)
+ real(kind=8),dimension(3) :: vectube
real(kind=8) :: Etube,xtemp,xminact,yminact,&
ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
sc_aa_tube,sc_bb_tube
!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
!C simple Kihara potential
subroutine calctube2(Etube)
- real(kind=8) :: vectube(3),enetube(nres*2)
+ real(kind=8),dimension(3) :: vectube
real(kind=8) :: Etube,xtemp,xminact,yminact,&
ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
sstube,ssgradtube,sc_aa_tube,sc_bb_tube
end subroutine calctube2
!=====================================================================================================================================
subroutine calcnano(Etube)
- real(kind=8) :: vectube(3),enetube(nres*2), &
- enecavtube(nres*2)
+ real(kind=8),dimension(3) :: vectube
+
real(kind=8) :: Etube,xtemp,xminact,yminact,&
ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact
integer:: i,j,iti
Etube=0.0d0
- print *,itube_start,itube_end,"poczatek"
+! print *,itube_start,itube_end,"poczatek"
do i=itube_start,itube_end
enetube(i)=0.0d0
enetube(i+nres)=0.0d0
!C fac=fac+faccav
!C 667 continue
endif
-
+ if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i)
do j=1,3
gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
enddo
+ if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres)
enddo
enddo
return
end subroutine set_shield_fac2
+!----------------------------------------------------------------------------
+! SOUBROUTINE FOR AFM
+ subroutine AFMvel(Eafmforce)
+ use MD_data, only:totTafm
+ real(kind=8),dimension(3) :: diffafm
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+!C Only for check grad COMMENT if not used for checkgrad
+!C totT=3.0d0
+!C--------------------------------------------------------
+!C print *,"wchodze"
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! totTafm=3.0
+ Eafmforce=0.5d0*forceAFMconst &
+ *(distafminit+totTafm*velAFMconst-afmdist)**2
+!C Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ enddo
+! print *,'AFM',Eafmforce,totTafm*velAFMconst,afmdist
+ return
+ end subroutine AFMvel
+!---------------------------------------------------------
+ subroutine AFMforce(Eafmforce)
+
+ real(kind=8),dimension(3) :: diffafm
+! real(kind=8) ::afmdist
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! print *,afmdist,distafminit
+ Eafmforce=-forceAFMconst*(afmdist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/afmdist
+ enddo
+!C print *,'AFM',Eafmforce
+ return
+ end subroutine AFMforce
!-----------------------------------------------------------------------------
#ifdef WHAM
allocate(grad_shield(3,-1:nres))
allocate(gg_tube_sc(3,-1:nres))
allocate(gg_tube(3,-1:nres))
+ allocate(gradafm(3,-1:nres))
!(3,maxres)
allocate(grad_shield_side(3,50,nres))
allocate(grad_shield_loc(3,50,nres))
! enddo
! enddo
- if (nss.gt.0) then
- allocate(idssb(nss),jdssb(nss))
+! if (nss.gt.0) then
+ allocate(idssb(maxdim),jdssb(maxdim))
+! allocate(newihpb(nss),newjhpb(nss))
!(maxdim)
- endif
+! endif
allocate(ishield_list(nres))
allocate(shield_list(50,nres))
allocate(dyn_ss_mask(nres))
allocate(fac_shield(nres))
+ allocate(enetube(nres*2))
+ allocate(enecavtube(nres*2))
+
!(maxres)
dyn_ss_mask(:)=.false.
!----------------------
projects / unres4.git / blobdiff