From 8d2ab9ba185dbbc31bfe3d1c66d7e1c9d632463b Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 31 Jul 2017 18:18:38 +0200 Subject: [PATCH] triss; AFM; Lorentz restrains included -debug might be on --- source/unres/MD.f90 | 6 + source/unres/control.F90 | 10 + source/unres/data/MD_data.f90 | 2 +- source/unres/data/control_data.f90 | 4 +- source/unres/data/energy_data.f90 | 19 +- source/unres/energy.f90 | 567 ++++++++++++++++++++++++++++++++---- source/unres/io.f90 | 113 ++++++- source/unres/io_base.f90 | 123 +++++++- source/unres/io_config.f90 | 7 + source/unres/unres.f90 | 3 +- 10 files changed, 761 insertions(+), 93 deletions(-) diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index 4368ce2..3725902 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -548,6 +548,7 @@ potE=potEcomp(0)-potEcomp(20) call cartgrad totT=totT+d_time + totTafm=totT ! Calculate the kinetic and total energy and the kinetic temperature call kinetic(EK) #ifdef MPI @@ -1001,6 +1002,7 @@ if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE,t_bath + totTafm=totT ! AL 4/17/17: Now writing d_t(0,:) too do i=0,2*nres write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) @@ -1287,6 +1289,7 @@ endif if (rattle) call rattle2 totT=totT+d_time + totTafm=totT if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 @@ -1726,6 +1729,7 @@ potE=potEcomp(0)-potEcomp(20) ! potE=energia_short(0)+energia_long(0) totT=totT+d_time + totTafm=totT ! Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -2403,6 +2407,7 @@ endif call random_vel totT=0.0d0 + totTafm=totT endif else ! Generate initial velocities @@ -2410,6 +2415,7 @@ write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 + totTafm=totT endif ! rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then diff --git a/source/unres/control.F90 b/source/unres/control.F90 index f45a88e..cc50b02 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -809,6 +809,14 @@ else call int_bounds(ndih_constr,idihconstr_start,idihconstr_end) endif + if (ntheta_constr.eq.0) then + ithetaconstr_start=1 + ithetaconstr_end=0 + else + call int_bounds & + (ntheta_constr,ithetaconstr_start,ithetaconstr_end) + endif + ! nsumgrad=(nres-nnt)*(nres-nnt+1)/2 ! nlen=nres-nnt+1 nsumgrad=(nres-nnt)*(nres-nnt+1)/2 @@ -1314,6 +1322,8 @@ iphi1_end=nres idihconstr_start=1 idihconstr_end=ndih_constr + ithetaconstr_start=1 + ithetaconstr_end=ntheta_constr iphid_start=iphi_start iphid_end=iphi_end-1 itau_start=4 diff --git a/source/unres/data/MD_data.f90 b/source/unres/data/MD_data.f90 index 73246ba..2bfe6e5 100644 --- a/source/unres/data/MD_data.f90 +++ b/source/unres/data/MD_data.f90 @@ -49,7 +49,7 @@ ntwx,ntwe logical :: mdpdb,large,print_compon,tbf,rest ! common /MDcalc/ - real(kind=8) :: totT,totE,potE,EK,amax,edriftmax,kinetic_T + real(kind=8) :: totT,totE,potE,EK,amax,edriftmax,kinetic_T,totTafm real(kind=8),dimension(:),allocatable :: potEcomp !(0:n_ene+4) ! common /lagrange/ real(kind=8),dimension(:,:),allocatable :: d_t,d_a,d_t_old !(3,0:MAXRES2) diff --git a/source/unres/data/control_data.f90 b/source/unres/data/control_data.f90 index 1b4e429..5f7d47e 100644 --- a/source/unres/data/control_data.f90 +++ b/source/unres/data/control_data.f90 @@ -30,11 +30,11 @@ ! common /cntrl/ integer :: modecalc,iscode,indpdb,indback,indphi,iranconf,& icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr,shield_mode,& - tubemode + tubemode,genconstr,afmlog,selfguide logical :: minim,refstr,pdbref,overlapsc,& energy_dec,sideadd,lsecondary,read_cart,unres_pdb,& vdisulf,searchsc,lmuca,dccart,extconf,out1file,& - gnorm_check,gradout,split_ene + gnorm_check,gradout,split_ene,with_theta_constr #ifdef CLUSTER integer :: iopt,nend,nstart,outpdb,outmol2 !cluster logical :: punch_dist,print_dist,lside,lprint_cart,lprint_int,& diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index 295d42a..71015c7 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -179,7 +179,7 @@ integer :: ns,nss,nfree integer,dimension(:),allocatable :: iss !(maxss) ! common /links/ - real(kind=8),dimension(:),allocatable :: dhpb,forcon,dhpb1 !(maxdim) !el dhpb1 !!! nie używane + real(kind=8),dimension(:),allocatable :: dhpb,forcon,dhpb1,fordepth !(maxdim) !el dhpb1 !!! nie używane integer :: nhpb integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane ! common /restraints/ @@ -187,7 +187,7 @@ ! common /links_split/ integer :: link_start,link_end ! common /dyn_ssbond/ - real(kind=8) :: Ht + real(kind=8) :: Ht,atriss,btriss,ctriss,dtriss integer,dimension(:),allocatable :: idssb,jdssb !(maxdim) logical :: dyn_ss logical,dimension(:),allocatable :: dyn_ss_mask !(maxres) @@ -212,12 +212,15 @@ !----------------------------------------------------------------------------- ! common.torcnstr ! common /torcnstr/ - integer :: ndih_constr,ndih_nconstr - integer,dimension(:),allocatable :: idih_constr,idih_nconstr !(maxdih_constr) - integer :: idihconstr_start,idihconstr_end + integer :: ndih_constr,ndih_nconstr,ntheta_constr + integer,dimension(:),allocatable :: idih_constr,idih_nconstr,itheta_constr !(maxdih_constr) + integer :: idihconstr_start,idihconstr_end, & + ithetaconstr_start,ithetaconstr_end real(kind=8) :: ftors - real(kind=8),dimension(:),allocatable :: drange !(maxdih_constr) + real(kind=8),dimension(:),allocatable :: drange,theta_constr0,theta_drange !(maxdih_constr) real(kind=8),dimension(:),allocatable :: phi0 !(maxdih_constr) + real(kind=8),dimension(:),allocatable :: for_thet_constr !(maxdih_constr) + !----------------------------------------------------------------------------- ! common.torsion ! Torsional constants of the rotation about virtual-bond dihedral angles @@ -264,7 +267,9 @@ real(kind=8),dimension(:), allocatable :: long_r_sidechain, & short_r_sidechain real(kind=8) :: VSolvSphere,VSolvSphere_div,buff_shield - +! AFM + real(kind=8) :: distafminit,forceAFMconst,velAFMconst + integer :: afmend,afmbeg diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index c166edf..565c695 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -123,7 +123,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 +222,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 +308,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 +371,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 +434,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 +521,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 +580,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 +624,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 +685,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 +694,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 +828,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 +859,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 +869,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 +892,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 +906,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 +928,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 +1401,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 +1421,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 +2795,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 +4898,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 !----------------------------------------------------------------------------- @@ -5328,7 +5427,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 +5453,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 +5630,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 +10274,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 +10300,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 +10453,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 +10486,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 +12842,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 +13132,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 +15436,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 +15591,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 +15677,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 +16095,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 +17847,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 +18163,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 +18185,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 +18221,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))& @@ -18348,7 +18733,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -18735,6 +19120,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 +19423,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,10 +19543,11 @@ 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)) diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 40bcb18..379978e 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -487,7 +487,7 @@ subroutine statout(itime) use energy_data - use control_data, only:refstr + use control_data use MD_data use MPI_data use compare, only:rms_nac_nnc @@ -513,8 +513,8 @@ character(len=4) :: format1,format2 character(len=30) :: format !el local variables - integer :: i,ii1,ii2 - real(kind=8) :: rms,frac,frac_nn,co + integer :: i,ii1,ii2,j + real(kind=8) :: rms,frac,frac_nn,co,distance #ifdef AIX if(itime.eq.0) then @@ -527,6 +527,46 @@ open(istat,file=statname,access="append") #endif #endif + if (AFMlog.gt.0) then + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')& + itime,totT,EK,potE,totE,& + rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& + potEcomp(23),me + format1="a133" + else +!C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,7f12.3,i5,$)') & + itime,totT,EK,potE,totE,& + kinetic_T,t_bath,gyrate(),& + potEcomp(23),me + format1="a114" + endif + else if (selfguide.gt.0) then + distance=0.0 + do j=1,3 + distance=distance+(c(j,afmend)-c(j,afmbeg))**2 + enddo + distance=dsqrt(distance) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, & + f9.3,i5,$)') & + itime,totT,EK,potE,totE,& + rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& + distance,potEcomp(23),me + format1="a133" +!C print *,"CHUJOWO" + else +!C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,8f12.3,i5,$)')& + itime,totT,EK,potE,totE, & + kinetic_T,t_bath,gyrate(),& + distance,potEcomp(23),me + format1="a114" + endif + else if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') & @@ -539,6 +579,7 @@ amax,kinetic_T,t_bath,gyrate(),me format1="a114" endif + ENDIF if(usampl.and.totT.gt.eq_time) then write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,& (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),& @@ -695,7 +736,7 @@ integer,dimension(maxres) :: itype_alloc integer :: iti,nsi,maxsi,itrial,itmp - real(kind=8) :: wlong,scalscp,co + real(kind=8) :: wlong,scalscp,co,ssscale allocate(weights(n_ene)) !----------------------------- allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres @@ -845,18 +886,23 @@ call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"ATRISS",atriss,0.301D0) + call reada(weightcard,"BTRISS",btriss,0.021D0) + call reada(weightcard,"CTRISS",ctriss,1.001D0) + call reada(weightcard,"DTRISS",dtriss,1.001D0) + call reada(weightcard,"SSSCALE",ssscale,1.0D0) dyn_ss=(index(weightcard,'DYN_SS').gt.0) call reada(weightcard,"HT",Ht,0.0D0) if (dyn_ss) then - ss_depth=ebr/wsc-0.25*eps(1,1) - Ht=Ht/wsc-0.25*eps(1,1) - akcm=akcm*wstrain/wsc - akth=akth*wstrain/wsc - akct=akct*wstrain/wsc - v1ss=v1ss*wstrain/wsc - v2ss=v2ss*wstrain/wsc - v3ss=v3ss*wstrain/wsc + ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale + Ht=(Ht/wsc-0.25*eps(1,1))*ssscale + akcm=akcm/wsc*ssscale + akth=akth/wsc*ssscale + akct=akct/wsc*ssscale + v1ss=v1ss/wsc*ssscale + v2ss=v2ss/wsc*ssscale + v3ss=v3ss/wsc*ssscale else ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain endif @@ -1023,6 +1069,7 @@ if (ndih_constr.gt.0) then allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr) allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr) + read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) if(me.eq.king.or..not.out1file)then @@ -1044,6 +1091,47 @@ phibound(2,ii) = phi0(i)+drange(i) enddo endif + if (with_theta_constr) then +!C with_theta_constr is keyword allowing for occurance of theta constrains + read (inp,*) ntheta_constr +!C ntheta_constr is the number of theta constrains + if (ntheta_constr.gt.0) then +!C read (inp,*) ftors + allocate(itheta_constr(ntheta_constr)) + allocate(theta_constr0(ntheta_constr)) + allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr)) + read (inp,*) (itheta_constr(i),theta_constr0(i), & + theta_drange(i),for_thet_constr(i), & + i=1,ntheta_constr) +!C the above code reads from 1 to ntheta_constr +!C itheta_constr(i) residue i for which is theta_constr +!C theta_constr0 the global minimum value +!C theta_drange is range for which there is no energy penalty +!C for_thet_constr is the force constant for quartic energy penalty +!C E=k*x**4 + if(me.eq.king.or..not.out1file)then + write (iout,*) & + 'There are',ntheta_constr,' constraints on phi angles.' + do i=1,ntheta_constr + write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), & + theta_drange(i), & + for_thet_constr(i) + enddo + endif + do i=1,ntheta_constr + theta_constr0(i)=deg2rad*theta_constr0(i) + theta_drange(i)=deg2rad*theta_drange(i) + enddo +!C if(me.eq.king.or..not.out1file) +!C & write (iout,*) 'FTORS',ftors +!C do i=1,ntheta_constr +!C ii = itheta_constr(i) +!C thetabound(1,ii) = phi0(i)-drange(i) +!C thetabound(2,ii) = phi0(i)+drange(i) +!C enddo + endif ! ntheta_constr.gt.0 + endif! with_theta_constr + nnt=1 #ifdef MPI if (me.eq.king) then @@ -1134,6 +1222,7 @@ call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb + if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co diff --git a/source/unres/io_base.f90 b/source/unres/io_base.f90 index 0a9dc14..0e1a986 100644 --- a/source/unres/io_base.f90 +++ b/source/unres/io_base.f90 @@ -68,6 +68,7 @@ write(iout,*) ' iss:',(iss(i),i=1,ns) ! Check whether the specified bridging residues are cystines. do i=1,ns + write(iout,*) i,iss(i) if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ',& @@ -83,6 +84,10 @@ #endif endif enddo + if (dyn_ss) then + if(.not.allocated(ihpb)) allocate(ihpb(ns)) + if(.not.allocated(jhpb)) allocate(jhpb(ns)) + endif ! Read preformed bridges. if (ns.gt.0) then read (inp,*) nss @@ -419,6 +424,30 @@ card=card(:ilen(card)+1)//karta return end subroutine card_concat +!---------------------------------------------------------------------------- + subroutine read_afminp + use geometry_data + use energy_data + use control_data, only:out1file + use MPI_data + character*320 afmcard + integer i + print *, "wchodze" + call card_concat(afmcard,.true.) + call readi(afmcard,"BEG",afmbeg,0) + call readi(afmcard,"END",afmend,0) + call reada(afmcard,"FORCE",forceAFMconst,0.0d0) + call reada(afmcard,"VEL",velAFMconst,0.0d0) + print *,'FORCE=' ,forceAFMconst +!------ NOW PROPERTIES FOR AFM + distafminit=0.0d0 + do i=1,3 + distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit + enddo + distafminit=dsqrt(distafminit) + print *,'initdist',distafminit + return + end subroutine read_afminp !----------------------------------------------------------------------------- subroutine read_dist_constr use MPI_data @@ -449,6 +478,16 @@ ! write (iout,*) "Calling read_dist_constr" ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup ! call flush(iout) + if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) + if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) + if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) + if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim)) + if(.not.allocated(forcon)) allocate(forcon(maxdim)) + if(.not.allocated(fordepth)) allocate(fordepth(maxdim)) + if ((genconstr.gt.0).and.(constr_dist.eq.11)) then + call gen_dist_constr2 + go to 1712 + endif call card_concat(controlcard,.true.) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) @@ -467,11 +506,11 @@ ! do i=1,npair_ ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) ! enddo - if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) - if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) - if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) - if(.not.allocated(forcon)) allocate(forcon(maxdim)) - +! if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) +! if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) +! if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) +! if(.not.allocated(forcon)) allocate(forcon(maxdim)) + call flush(iout) do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup @@ -546,10 +585,29 @@ endif enddo do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +! if (forcon(nhpb+1).gt.0.0d0) then +! nhpb=nhpb+1 +! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & + ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + else +!C print *,"in else" + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & + ibecarb(i),forcon(nhpb+1) + endif if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) & + dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif + #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& @@ -558,12 +616,60 @@ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif - endif enddo + 1712 continue call flush(iout) return end subroutine read_dist_constr !----------------------------------------------------------------------------- + subroutine gen_dist_constr2 + use MPI_data + ! use control + use geometry, only: dist + use geometry_data + use control_data + use energy_data + integer :: i,j + real(kind=8) :: distance + if (constr_dist.eq.11) then + do i=nstart_sup,nstart_sup+nsup-1 + do j=i+2,nstart_sup+nsup-1 + distance=dist(i,j) + if (distance.le.15.0) then + nhpb=nhpb+1 + ihpb(nhpb)=i+nstart_seq-nstart_sup + jhpb(nhpb)=j+nstart_seq-nstart_sup + forcon(nhpb)=sqrt(0.04*distance) + fordepth(nhpb)=sqrt(40.0/distance) + dhpb(nhpb)=distance-0.1d0 + dhpb1(nhpb)=distance+0.1d0 + +#ifdef MPI + if (.not.out1file .or. me.eq.king) & + write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", & + nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#else + write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", & + nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#endif + endif + enddo + enddo + else + do i=nstart_sup,nstart_sup+nsup-1 + do j=i+2,nstart_sup+nsup-1 + nhpb=nhpb+1 + ihpb(nhpb)=i+nstart_seq-nstart_sup + jhpb(nhpb)=j+nstart_seq-nstart_sup + forcon(nhpb)=weidis + dhpb(nhpb)=dist(i,j) + enddo + enddo + endif + return + end subroutine gen_dist_constr2 + +!----------------------------------------------------------------------------- #ifdef WINIFL subroutine flush(iu) return @@ -1130,6 +1236,7 @@ ! iii=igeom igeom=iout #endif + print *,nss IF (NSS.LE.9) THEN #ifdef CLUSTER WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS) diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 91daac1..d3c0368 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -2970,6 +2970,12 @@ !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "shield_mode",shield_mode !C Varibles set size of box + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr + AFMlog=(index(controlcard,'AFM')) + selfguide=(index(controlcard,'SELFGUIDE')) + print *,'AFMlog',AFMlog,selfguide,"KUPA" + call readi(controlcard,'GENCONSTR',genconstr,0) call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) @@ -4052,6 +4058,7 @@ open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath + totTafm=totT ! do i=1,2*nres ! AL 4/17/17: Now reading d_t(0,:) too do i=0,2*nres diff --git a/source/unres/unres.f90 b/source/unres/unres.f90 index 619aa73..40c0140 100644 --- a/source/unres/unres.f90 +++ b/source/unres/unres.f90 @@ -348,11 +348,12 @@ else print *,'refstr=',refstr if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + print *,"after rms_nac_ncc" call briefout(0,etot) endif if (outpdb) call pdbout(etot,titel(:32),ipdb) if (outmol2) call mol2out(etot,titel(:32)) -!elwrite(iout,*) "after exec_eeval_or_minim" + write(iout,*) "after exec_eeval_or_minim" return end subroutine exec_eeval_or_minim !----------------------------------------------------------------------------- -- 1.7.9.5