PBC working for EGB (only)
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 25 Apr 2017 12:06:30 +0000 (14:06 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 25 Apr 2017 12:06:30 +0000 (14:06 +0200)
12 files changed:
.gitignore
PARAM/bond_AM1_ext_dum.parm [new file with mode: 0644]
source/unres/MCM_MD.f90
source/unres/data/control_data.f90
source/unres/data/energy_data.f90
source/unres/data/geometry_data.f90
source/unres/energy.f90
source/unres/geometry.f90
source/unres/io_config.f90
source/unres/minim.f90
source/unres/regularize.f90
source/unres/unres.f90

index f906e40..99fb1e8 100644 (file)
@@ -16,6 +16,7 @@ bin/*exe
 build/
 build2/
 build*/
+run/
 
 # latex files in documentation 
 doc/*/*.aux
diff --git a/PARAM/bond_AM1_ext_dum.parm b/PARAM/bond_AM1_ext_dum.parm
new file mode 100644 (file)
index 0000000..cf2100e
--- /dev/null
@@ -0,0 +1,26 @@
+1 3.800 1.900  41.7  0.000                                  43.0  14.33  2.5  ! peptide group and dummy atom
+1 1.396 243.1  0.000                                        59.0  19.67  5.0  ! Cys
+2 2.103  71.3  0.850 2.500 128.0  0.033                     88.0  29.33  6.2  ! Met
+1 2.997 124.6  0.000                                       104.0  34.67  6.8  ! Phe
+2 1.645 260.7  0.857 1.908 312.6 .00010                     70.0  23.33  6.2  ! Ile
+2 1.782 638.3  2.360 2.086 160.4  0.409                     70.0  23.33  6.3  ! Leu
+1 1.488 294.4  0.000                                        56.0  18.67  5.8  ! Val
+2 3.368 123.1  0.000 3.686 129.1 .00049                    143.0  47.67  7.2  ! Trp
+1 3.362 113.2  0.000                                       120.0  40.00  6.9  ! Tyr
+1 0.778 353.0  0.000                                        28.0   9.33  4.6  ! Ala
+0                                                           14.0   0.00  3.8  ! Gly
+1 1.480 295.8  0.000                                        58.0  19.33  5.6  ! Thr
+1 1.311 269.6  0.000                                        44.0  14.67  4.8  ! Ser
+3 2.125 147.9  2.125 2.424 138.9  .0333 2.776 383.8  0.000  85.0  28.33  6.1  ! Gln
+1 2.008 161.4  0.000                                        71.0  23.67  5.7  ! Asn
+3 2.093 131.2  1.943 2.425 146.5  .0263 2.784 479.3  0.784  85.0  28.33  6.1  ! Glu
+1 2.030 160.2  0.000                                        71.0  23.67  5.6  ! Asp
+1 2.739 134.5  0.000                                        95.0  31.67  6.2  ! His
+3 2.644  48.8  1.707 3.433  34.8 .00123 4.080 899.9  1.175 114.0  38.00  6.8  ! Arg
+3 2.379  99.0  1.974 2.704 157.5  0.546 3.073 164.7  0.055  86.0  28.67  6.3  ! Lys
+1 1.422 605.2  0.000                                        71.0  23.67  5.6  ! Pro
+2 2.142  71.3  0.850 2.500 128.0  0.033                    135.0  45.00  6.2  ! SeMet
+1 3.799 124.6  0.000                                       149.0  49.67  7.2  ! Dap(Bz)
+1 0.743 353.00 0.000                                        42.0  14.00  4.7  ! Aib
+1 1.210 353.00 0.000                                        42.0  14.00  5.6  ! Abu
+
index aa3f3fe..c518414 100644 (file)
       nlist=0
 #ifdef UNRES
       call var_to_geom(nvar,x)
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
       do k=1,2*nres
        do kk=1,3
 
 #ifdef UNRES
        call var_to_geom(nvar,x1)
+      write(iout,*) 'Warning calling chainbuild'
        call chainbuild
 !d     write(iout,*)'C and CREF'
 !d     write(iout,'(i5,3f10.5,5x,3f10.5)')(k,(c(j,k),j=1,3),
       WhatsUp=0
       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
       write (iout,'(/80(1h*)/a)') 'Initial energies:'
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
       call etotal(energia)
       etot = energia(0)
 !d        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
           call var_to_geom(nvar,varia)
 ! Rebuild the chain.
+          write(iout,*) 'Warning calling chainbuild'
           call chainbuild
           MoveType=0
           nbond=0
                  MoveType,' returned from PERTURB.'
               goto 20
             endif
+            write(iout,*) 'Warning calling chainbuild'
             call chainbuild
           else
             MoveType=0
       WhatsUp=0
       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
       write (iout,'(/80(1h*)/a)') 'Initial energies:'
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
       call geom_to_var(nvar,varia)
       call etotal(energia)
           enddo
           call var_to_geom(nvar,varia)
 ! Rebuild the chain.
+          write(iout,*) 'Warning calling chainbuild'
           call chainbuild
           MoveType=0
           nbond=0
               varia(i)=xpool(i,ii)
             enddo
             call var_to_geom(nvar,varia)
+            write(iout,*) 'Warning calling chainbuild'
             call chainbuild  
 !d          call intout
 !d          write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
index 1e2bbf9..414ac0c 100644 (file)
@@ -50,7 +50,7 @@
 !-----------------------------------------------------------------------------
 ! common.spitele
 !      common /splitele/
-      real(kind=8) :: r_cut,rlamb
+      real(kind=8) :: r_cut,rlamb,r_cut_ele,rlamb_ele
 !-----------------------------------------------------------------------------
 ! common.time1
 !     FOUND_NAN - set by calcf to stop sumsl via stopx
index 873f053..1b05b85 100644 (file)
       real(kind=8),dimension(:,:),allocatable :: r0d,eps_scp,rscp !(ntyp,2)  r0d  !!! nie używane
 ! 12/5/03 modified 09/18/03 Bond stretching parameters.
 !      common /stretch/
-      real(kind=8) :: vbldp0,akp,distchainmax
+      real(kind=8) :: vbldp0,akp,distchainmax,vbldpDUM
       real(kind=8),dimension(:,:),allocatable :: vbldsc0,aksc,abond0 !(maxbondterm,ntyp)
       integer,dimension(:),allocatable :: nbondterm    !(ntyp)
 !-----------------------------------------------------------------------------
 !      common/fourier/  z wham
       real(kind=8),dimension(:,:),allocatable :: b !(13,0:maxtor)
 !-----------------------------------------------------------------------------
+! 24 Apr 2017 
+! Varibles for cutoff on electorstatic
+      real(kind=8) sss_ele_cut,sss_ele_grad
+      integer xshift,yshift,zshift
 !-----------------------------------------------------------------------------
+
       end module energy_data
index 5efd1be..022d4ab 100644 (file)
@@ -81,4 +81,7 @@
       integer,dimension(:),allocatable :: itype_pdb !(maxres) initialize in molread
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
+!  common.box
+      real(kind=8) :: boxxsize,boxysize,boxzsize
+
       end module geometry_data
index 40153e4..1d58136 100644 (file)
 
 #ifdef MPI      
       real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
+! shielding effect varibles for MPI
+!      real(kind=8)   fac_shieldbuf(maxres),
+!     & grad_shield_locbuf(3,maxcontsshi,-1:maxres),
+!     & grad_shield_sidebuf(3,maxcontsshi,-1:maxres),
+!     & grad_shieldbuf(3,-1:maxres)
+!       integer ishield_listbuf(maxres),
+!     &shield_listbuf(maxcontsshi,maxres)
+
 !      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 !     & " nfgtasks",nfgtasks
       if (nfgtasks.gt.1) then
 !      include 'COMMON.SBRIDGE'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj
+      integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
       real(kind=8) :: evdw,sig0ij
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
       integer :: ii
 !cccc      energy_dec=.false.
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          xi=dmod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=dmod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=dmod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
 !           alf1=0.0D0
 !           alf2=0.0D0
 !           alf12=0.0D0
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+           xj=c(1,nres+j)
+           yj=c(2,nres+j)
+           zj=c(3,nres+j)
+          xj=dmod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=dmod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=dmod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=0
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+       else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+       endif
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
 !          write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+!            print *,sss_ele_cut,sss_ele_grad,&
+!            1.0d0/(rij),r_cut_ele,rlamb_ele
+            if (sss_ele_cut.le.0.0) cycle
 ! Calculate angle-dependent terms of energy and contributions to their
 ! derivatives.
             call sc_angular
 !          write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d
 !          " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d
             evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij
+            evdw=evdw+evdwij*sss_ele_cut
             if (lprn) then
             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac
+!            print *,'before fac',fac,rij,evdwij
+            fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij
+!            print *,'grad part scale',fac,   &
+!             evdwij*sss_ele_grad/sss_ele_cut &
+!            /sigma(itypi,itypj)*rij
 !            fac=0.0d0
 ! Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
+!            print *,'before sc_grad', gg(1),gg(2),gg(3)
 ! Calculate angular part of the gradient.
             call sc_grad
             ENDIF    ! dyn_ss            
 !      if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
 
       do i=ibondp_start,ibondp_end
+        if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
-          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
-          do j=1,3
-          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
-            *dc(j,i-1)/vbld(i)
-          enddo
-          if (energy_dec) write(iout,*) &
-             "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+!C          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+!C          do j=1,3
+!C          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
+!C            *dc(j,i-1)/vbld(i)
+!C          enddo
+!C          if (energy_dec) write(iout,*) &
+!C             "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+        diff = vbld(i)-vbldpDUM
         else
         diff = vbld(i)-vbldp0
+        endif
         if (energy_dec) write (iout,'(a7,i5,4f7.3)') &
            "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
         enddo
 !        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
-        endif
+!        endif
       enddo
       estr=0.5d0*AKP*estr+estr1
 !
 !       " sigder",sigder
 !      write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
 !      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+!C      print *,sss_ele_cut,'in sc_grad'
       do k=1,3
         dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
       enddo
       do k=1,3
-        gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
+!C      print *,'gg',k,gg(k)
       enddo 
 !      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
         gvdwx(k,i)=gvdwx(k,i)-gg(k) &
                   +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
-                  +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+                  +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv    &
+                  *sss_ele_cut
+
         gvdwx(k,j)=gvdwx(k,j)+gg(k) &
                   +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
-                  +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+                  +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv    &
+                  *sss_ele_cut
+
 !        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
 !                 +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 !        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
@@ -10519,7 +10594,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !      call intcartderiv
 !      call checkintcartgrad
       call zerograd
-      aincr=1.0D-4
+      aincr=1.0D-6
       write(iout,*) 'Calling CHECK_ECARTINT.'
       nf=0
       icall=0
@@ -10924,6 +10999,35 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       endif
       return
       end function sscale
+!!!!!!!!!! PBCSCALE
+      real(kind=8) function sscale_ele(r)
+!      include "COMMON.SPLITELE"
+      real(kind=8) :: r,gamm
+      if(r.lt.r_cut_ele-rlamb_ele) then
+        sscale_ele=1.0d0
+      else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+        gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+        sscale_ele=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+      else
+        sscale_ele=0d0
+      endif
+      return
+      end function sscale_ele
+
+      real(kind=8)  function sscagrad_ele(r)
+      real(kind=8) :: r,gamm
+!      include "COMMON.SPLITELE"
+      if(r.lt.r_cut_ele-rlamb_ele) then
+        sscagrad_ele=0.0d0
+      else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+        gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+        sscagrad_ele=gamm*(6*gamm-6.0d0)/rlamb_ele
+      else
+        sscagrad_ele=0.0d0
+      endif
+      return
+      end function sscagrad_ele
+!!!!!!!!!!!!!!!
 !-----------------------------------------------------------------------------
       subroutine elj_long(evdw)
 !
@@ -11542,9 +11646,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !      include 'COMMON.CONTROL'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj
+      integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
       real(kind=8) :: sss,e1,e2,evdw
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
+
       evdw=0.0D0
 !cccc      energy_dec=.false.
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -11559,6 +11666,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -11589,16 +11702,57 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+! Searching for nearest neighbour
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          xj_safe=xj
+          yj_safe=yj
+          zj_safe=zj
+          subchap=0
+          do xshift=-1,1
+          do yshift=-1,1
+          do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+          enddo
+          enddo
+          enddo
+          if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+          else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+          endif
+
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            if (sss_ele_cut.le.0.0) cycle
             if (sss.lt.1.0d0) then
 
 ! Calculate angle-dependent terms of energy and contributions to their
@@ -11629,7 +11783,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 !     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
-              evdw=evdw+evdwij*(1.0d0-sss)
+              evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut
               if (lprn) then
               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -11651,6 +11805,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
+              fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -11687,9 +11843,11 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !      include 'COMMON.CONTROL'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj
+      integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
       real(kind=8) :: sss,e1,e2,evdw,rij_shift
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
       evdw=0.0D0
 !cccc      energy_dec=.false.
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -11704,6 +11862,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+        dxi=dc_norm(1,nres+i)
+        dyi=dc_norm(2,nres+i)
+        dzi=dc_norm(3,nres+i)
+!        dsci_inv=dsc_inv(itypi)
+        dsci_inv=vbld_inv(i+nres)
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -11734,15 +11904,60 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+!            xj=c(1,nres+j)-xi
+!            yj=c(2,nres+j)-yi
+!            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+! Searching for nearest neighbour
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          xj_safe=xj
+          yj_safe=yj
+          zj_safe=zj
+          subchap=0
+          do xshift=-1,1
+          do yshift=-1,1
+          do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+          enddo
+          enddo
+          enddo
+          if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+          else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+          endif
+
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+            sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            if (sss_ele_cut.le.0.0) cycle
 
             if (sss.gt.0.0d0) then
 
@@ -11796,6 +12011,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
+              fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij
+
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -13291,16 +13509,19 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
       enddo
       do k=1,3
-        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac&
+         *sss_ele_cut
       enddo 
 !      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
         gvdwx(k,i)=gvdwx(k,i)-gg(k) &
                   +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
-                +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
+                +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac&
+                 *sss_ele_cut
         gvdwx(k,j)=gvdwx(k,j)+gg(k) &
                   +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
-                +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
+                +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac&
+         *sss_ele_cut
 !        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
 !     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 !        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
index 0f3febc..f9751f6 100644 (file)
       do i=1,nres-1
 !in wham      do i=1,nres
         iti=itype(i)
-        if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+        if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
+       (iti.ne.ntyp1  .and. itype(i+1).ne.ntyp1)) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 !test          stop
         endif
index 018dbee..da0414d 100644 (file)
         endif
       enddo
 #else
-      read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
+      read (ibond,*) junk,vbldpDUM,vbldp0,akp,rjunk,mp,ip,pstok
       do i=1,ntyp
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
          j=1,nbondterm(i)),msc(i),isc(i),restok(i)
           goto 10
         else if (card(:3).eq.'TER') then
 ! End current chain
-          ires_old=ires+1
+          ires_old=ires+2
           ishift1=ishift1+1
           itype(ires_old)=ntyp1
+          itype(ires_old-1)=ntyp1
           ibeg=2
 !          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
       nres=ires
       do i=2,nres-1
 !        write (iout,*) i,itype(i)
-        if (itype(i).eq.ntyp1) then
+!        if (itype(i).eq.ntyp1) then
 !          write (iout,*) "dummy",i,itype(i)
-          do j=1,3
-            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
+!          do j=1,3
+!            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
-            dc(j,i)=c(j,i)
-          enddo
-        endif
+!            dc(j,i)=c(j,i)
+!          enddo
+!        endif
+        if (itype(i).eq.ntyp1) then
+         if (itype(i+1).eq.ntyp1) then
+! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+           if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the last dummy residue
+!            print *,i,'tu dochodze'
+            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+            if (fail) then
+              e2(1)=0.0d0
+              e2(2)=1.0d0
+              e2(3)=0.0d0
+            endif !fail
+            print *,i,'a tu?'
+            do j=1,3
+             c(j,i)=c(j,i-1)-1.9d0*e2(j)
+            enddo
+           else   !unres_pdb
+           do j=1,3
+             dcj=(c(j,i-2)-c(j,i-3))/2.0
+            if (dcj.eq.0) dcj=1.23591524223
+             c(j,i)=c(j,i-1)+dcj
+             c(j,nres+i)=c(j,i)
+           enddo
+          endif   !unres_pdb
+         else     !itype(i+1).eq.ntyp1
+          if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+            if (fail) then
+              e2(1)=0.0d0
+              e2(2)=1.0d0
+              e2(3)=0.0d0
+            endif
+            do j=1,3
+              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+            enddo
+          else !unres_pdb
+           do j=1,3
+            dcj=(c(j,i+3)-c(j,i+2))/2.0
+            if (dcj.eq.0) dcj=1.23591524223
+            c(j,i)=c(j,i+1)-dcj
+            c(j,nres+i)=c(j,i)
+           enddo
+          endif !unres_pdb
+         endif !itype(i+1).eq.ntyp1
+        endif  !itype.eq.ntyp1
+
       enddo
 ! Calculate the CM of the last side chain.
       if (iii.gt.0)  then
       timem=timlim
       modecalc=0
       call reada(controlcard,"T_BATH",t_bath,300.0d0)
+!C  Varibles set size of box
+      call reada(controlcard,'BOXX',boxxsize,100.0d0)
+      call reada(controlcard,'BOXY',boxysize,100.0d0)
+      call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+! CUTOFFF ON ELECTROSTATICS
+      call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+      call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
+
+!C-------------------------
       minim=(index(controlcard,'MINIMIZE').gt.0)
       dccart=(index(controlcard,'CART').gt.0)
       overlapsc=(index(controlcard,'OVERLAP').gt.0)
index 4305640..d556363 100644 (file)
 !     else
 !       not_done=.false.
 !     endif
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
 
 !el---------------------
       endif
 #endif
       call var_to_geom_restr(n,x)
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild 
 !
 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
       icg=mod(nf,2)+1
       call var_to_geom_restr(n,x)
       call zerograd
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
 !d    write (iout,*) 'ETOTAL called from FUNC'
       call etotal(energia)
           cur_omeg(i)=omeg(i)
         endif
       enddo
-
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
       call egb1(evdw)
       call esc(escloc)
       wang=orig_w(11)
       wtor=orig_w(13)
       wtor_d=orig_w(14)
-
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
       call etotal(energy_)
       etot=energy_(0)
 
       call var_to_geom_restr(n,x)
       call zerograd
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild
 !d    write (iout,*) 'ETOTAL called from FUNC'
       call egb1(evdw)
       if (nf.eq.0) return
       goto 40
    30 call var_to_geom_restr(n,x)
+      write(iout,*) 'Warning calling chainbuild'
       call chainbuild 
 !
 ! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
index 0eb81e5..70214b7 100644 (file)
@@ -45,6 +45,7 @@
        ' nstart_seq=',nstart_seq,' nstart_sup',nstart_sup
       write (iout,'(/a/)') 'Initial energies:'
       call geom_to_var(nvar,varia)
+      write(iout,*) 'Warning: Calling chainbuild'
       call chainbuild
       call etotal(energia)
       etot=energia(0)
index 4216238..deb5713 100644 (file)
 
       integer :: j,k
       call alloc_compare_arrays
-      if (indpdb.eq.0) call chainbuild
+      if (indpdb.eq.0) then 
+      call chainbuild
+      write(iout,*) 'Warning: Calling chainbuild'
+      endif
 #ifdef MPI
       time00=MPI_Wtime()
 #endif
         else 
           if (indpdb.ne.0) then 
             call bond_regular
+            write(iout,*) 'Warning: Calling chainbuild'
             call chainbuild
           endif
           call geom_to_var(nvar,varia)
             read (intin,'(i5)',end=1100,err=1100) iconf
             call read_angles(intin,*11)
             call geom_to_var(nvar,varia)
+            write(iout,*) 'Warning: Calling chainbuild'
             call chainbuild
           endif
           write (iout,'(a,i7)') 'Conformation #',iconf
             read (intin,'(i5)',end=11,err=11) iconf
             call read_angles(intin,*11)
             call geom_to_var(nvar,varia)
+            write(iout,*) 'Warning: Calling chainbuild'
             call chainbuild
           endif
           write (iout,'(a,i7)') 'Conformation #',iconf
 !         print *,'result received from worker ',man,' sending now'
 
           call var_to_geom(nvar,varia)
+          write(iout,*) 'Warning: Calling chainbuild'
           call chainbuild
           call etotal(energy_)
           iconf=ind(2)
             read (intin,'(i5)',end=1101,err=1101) iconf
             call read_angles(intin,*11)
             call geom_to_var(nvar,varia)
+            write(iout,*) 'Warning: Calling chainbuild'
             call chainbuild
           endif
           n=n+1
                      CG_COMM,muster,ierr)
 
         call var_to_geom(nvar,varia)
+        write(iout,*) 'Warning: Calling chainbuild'
         call chainbuild
         call etotal(energy_)
         iconf=ind(2)
             read (intin,'(i5)',end=11,err=11) iconf
             call read_angles(intin,*11)
             call geom_to_var(nvar,varia)
+            write(iout,*) 'Warning: Calling chainbuild'
             call chainbuild
           endif
         write (iout,'(a,i7)') 'Conformation #',iconf
 !        if (itype(i).ne.10) 
 !     &      vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
 !      enddo
-      if (indpdb.eq.0) call chainbuild
+      if (indpdb.eq.0) then 
+        write(iout,*) 'Warning: Calling chainbuild'
+        call chainbuild
+      endif
 !      do i=0,nres
 !        do j=1,3
 !          dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
       logical :: debug
 
       call alloc_compare_arrays
+      write(iout,*) 'Warning: Calling chainbuild'
       call chainbuild
       call etotal(energy_)
       call enerprint(energy_)