X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=4065e96d235689cae964a25825b192caad81394a;hb=4fa6ed0fa1ee37552df6064fe60a73f218c101ca;hp=c166edfd1386976b8ebfee474f408f0a044792b9;hpb=affbf55b871a9bcc8fdb91efeb8e6bf7a3d8d003;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index c166edf..4065e96 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -45,6 +45,7 @@ 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/ @@ -123,7 +124,7 @@ 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) @@ -222,7 +223,8 @@ 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 @@ -307,7 +309,7 @@ #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. @@ -370,7 +372,7 @@ ! 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 & @@ -433,7 +435,7 @@ ! Calculate the virtual-bond-angle energy. ! if (wang.gt.0d0) then - call ebend(ebe) + call ebend(ebe,ethetacnstr) else ebe=0 endif @@ -520,6 +522,13 @@ 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 @@ -572,6 +581,8 @@ 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 @@ -614,7 +625,7 @@ 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 @@ -675,6 +686,8 @@ 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 & @@ -682,14 +695,17 @@ +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 @@ -813,7 +829,7 @@ 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) @@ -844,6 +860,8 @@ 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,& @@ -852,8 +870,8 @@ 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)'/ & @@ -875,9 +893,11 @@ '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 @@ -887,7 +907,8 @@ 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)'/ & @@ -908,9 +929,11 @@ '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 @@ -1379,7 +1402,7 @@ 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) @@ -1399,6 +1422,27 @@ '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)) @@ -2752,8 +2796,8 @@ ! 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 @@ -4855,45 +4899,101 @@ 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 !----------------------------------------------------------------------------- @@ -5085,6 +5185,9 @@ 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 @@ -5328,7 +5431,7 @@ 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. @@ -5354,7 +5457,10 @@ !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 @@ -5528,6 +5634,37 @@ 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 @@ -10141,6 +10278,7 @@ 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)& @@ -10166,6 +10304,7 @@ 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) & @@ -10318,6 +10457,7 @@ 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) & @@ -10350,6 +10490,7 @@ 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,) & @@ -12705,12 +12846,34 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -12973,6 +13136,28 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -15255,7 +15440,7 @@ write(iout,*) 'Calling CHECK_ECARTIN 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 @@ -15410,7 +15595,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! ! Calculate the virtual-bond-angle energy. ! - call ebend(ebe) + call ebend(ebe,ethetacnstr) ! ! Calculate the SC local energy. ! @@ -15496,7 +15681,35 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -15886,6 +16099,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -17637,6 +17851,176 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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] @@ -17783,15 +18167,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -17802,11 +18189,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -17836,7 +18225,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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))& @@ -17942,7 +18331,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 @@ -18103,7 +18492,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 @@ -18340,15 +18729,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -18441,7 +18830,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 @@ -18544,6 +18933,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -18735,6 +19125,62 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -18982,6 +19428,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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)) @@ -19101,14 +19548,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! 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. !----------------------