triss; AFM; Lorentz restrains included -debug might be on
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 31 Jul 2017 16:18:38 +0000 (18:18 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 31 Jul 2017 16:18:38 +0000 (18:18 +0200)
source/unres/MD.f90
source/unres/control.F90
source/unres/data/MD_data.f90
source/unres/data/control_data.f90
source/unres/data/energy_data.f90
source/unres/energy.f90
source/unres/io.f90
source/unres/io_base.f90
source/unres/io_config.f90
source/unres/unres.f90

index 4368ce2..3725902 100644 (file)
       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
         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)
           endif                    
           if (rattle) call rattle2
           totT=totT+d_time
+        totTafm=totT
           if (d_time.ne.d_time0) then
             d_time=d_time0
 #ifndef   LANG0
       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
         endif
         call random_vel
         totT=0.0d0
+        totTafm=totT
        endif
       else
 ! Generate initial velocities
          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
index f45a88e..cc50b02 100644 (file)
       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
       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
index 73246ba..2bfe6e5 100644 (file)
@@ -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)
index 1b4e429..5f7d47e 100644 (file)
 !      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,&
index 295d42a..71015c7 100644 (file)
       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/
 !      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)
 !-----------------------------------------------------------------------------
 ! 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
       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
 
 
 
index c166edf..565c695 100644 (file)
         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
 !-----------------------------------------------------------------------------
       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,) &
@@ -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))
index 40bcb18..379978e 100644 (file)
       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
       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
       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,$)') &
                  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),&
       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
       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
       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
           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
         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
index 0a9dc14..0e1a986 100644 (file)
@@ -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 ',&
 #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
       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
 !      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)
 !      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
         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 ",&
           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
 !      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)
index 91daac1..d3c0368 100644 (file)
 !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)
 
       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
index 619aa73..40c0140 100644 (file)
       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
 !-----------------------------------------------------------------------------