DO zmiany czytanie pisanie i cutoff

diff --git a/.gitignore b/.gitignore
index d5a5c53..77a69e0 100644 (file)
--- a/.gitignore
+++ b/.gitignore
@@ -27,3 +27,4 @@ bin/unres/MD/unres_ifort_MPICH_GAB_czyt.exe
 bin/unres/MD-M/unres_Tc_procor_newparm_em64-D-symetr.exe
 DIL/
 compinfo
+period_build

diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN
index f343887..bf5f688 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CHAIN
+++ b/source/unres/src_MD-M/COMMON.CHAIN
@@ -14,3 +14,5 @@
      & nsup,nstart_sup,nstart_seq,
      & chain_length,iprzes,tabperm(maxperm,maxsym),nperm
       common /from_zscore/ nz_start,nz_end,iz_sc
+      double precision boxxsize,boxysize,boxzsize
+      common /box/  boxxsize,boxysize,boxzsize 

diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT
index 83af3fb..37d3e88 100644 (file)
--- a/source/unres/src_MD-M/COMMON.INTERACT
+++ b/source/unres/src_MD-M/COMMON.INTERACT
@@ -24,8 +24,10 @@ C 12/1/95 Array EPS included in the COMMON block.
      & rpp(2,2),epp(2,2),elpp6(2,2),elpp3(2,2),eps_scp(ntyp,2),
      & rscp(ntyp,2)
 c 12/5/03 modified 09/18/03 Bond stretching parameters.
-      double precision vbldp0,vbldsc0,akp,aksc,abond0,distchainmax
+      double precision vbldp0,vbldpDUM,
+     &  vbldsc0,akp,aksc,abond0,distchainmax
       integer nbondterm
-      common /stretch/ vbldp0,vbldsc0(maxbondterm,ntyp),akp,
+      common /stretch/ vbldp0,vbldpDUM,
+     & vbldsc0(maxbondterm,ntyp),akp,
      & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp),
      & distchainmax,nbondterm(ntyp)

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index d4e2092..635a41e 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -1413,6 +1413,7 @@ C
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
       logical lprn
+      integer xshift,yshift,zshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -1420,6 +1421,11 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.false.
       ind=0
+C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0
+C we have the original box)
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
       do i=iatsc_s,iatsc_e
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
@@ -1427,6 +1433,32 @@ c     if (icall.eq.0) lprn=.false.
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+C Return atom into box, boxxsize is size of box in x dimension
+  134   continue
+        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
+        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+C Condition for being inside the proper box
+        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
+     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
+        go to 134
+        endif
+  135   continue
+        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
+        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+C Condition for being inside the proper box
+        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
+     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
+        go to 135
+        endif
+  136   continue
+        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
+        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+C Condition for being inside the proper box
+        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
+     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
+        go to 136
+        endif
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1467,12 +1499,41 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           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)
+C Return atom J into box the original box
+  137   continue
+        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+C Condition for being inside the proper box
+        if ((xj.gt.((0.5d0)*boxxsize)).or.
+     &       (xj.lt.((-0.5d0)*boxxsize))) then
+        go to 137
+        endif
+  138   continue
+        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+C Condition for being inside the proper box
+        if ((yj.gt.((0.5d0)*boxysize)).or.
+     &       (yj.lt.((-0.5d0)*boxysize))) then
+        go to 138
+        endif
+  139   continue
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+C Condition for being inside the proper box
+        if ((zj.gt.((0.5d0)*boxzsize)).or.
+     &       (zj.lt.((-0.5d0)*boxzsize))) then
+        go to 139
+        endif
+
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
+            xj=xj-xi
+            yj=yj-yi
+            zj=zj-zi
 c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
 c            write (iout,*) "j",j," dc_norm",
 c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
@@ -1536,6 +1597,9 @@ C Calculate angular part of the gradient.
           enddo      ! j
         enddo        ! iint
       enddo          ! i
+      enddo          ! zshift
+      enddo          ! yshift
+      enddo          ! xshift
 c      write (iout,*) "Number of loop steps in EGB:",ind
 cccc      energy_dec=.false.
       return
@@ -2788,6 +2852,7 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms
 C
 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
+C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
@@ -2800,6 +2865,31 @@ C
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+C Return atom into box, boxxsize is size of box in x dimension
+  184   continue
+        if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
+        if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+C Condition for being inside the proper box
+        if ((xmedi.gt.((0.5d0)*boxxsize)).or.
+     &       (xmedi.lt.((-0.5d0)*boxxsize))) then
+        go to 184
+        endif
+  185   continue
+        if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
+        if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+C Condition for being inside the proper box
+        if ((ymedi.gt.((0.5d0)*boxysize)).or.
+     &       (ymedi.lt.((-0.5d0)*boxysize))) then
+        go to 185
+        endif
+  186   continue
+        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
+        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+C Condition for being inside the proper box
+        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
+     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
+        go to 186
+        endif
         num_conti=0
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -2818,12 +2908,42 @@ C
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+C Return atom into box, boxxsize is size of box in x dimension
+  194   continue
+        if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
+        if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+C Condition for being inside the proper box
+        if ((xmedi.gt.((0.5d0)*boxxsize)).or.
+     &       (xmedi.lt.((-0.5d0)*boxxsize))) then
+        go to 194
+        endif
+  195   continue
+        if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
+        if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+C Condition for being inside the proper box
+        if ((ymedi.gt.((0.5d0)*boxysize)).or.
+     &       (ymedi.lt.((-0.5d0)*boxysize))) then
+        go to 195
+        endif
+  196   continue
+        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
+        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+C Condition for being inside the proper box
+        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
+     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
+        go to 196
+        endif
+
         num_conti=num_cont_hb(i)
         call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
      &   call eturn4(i,eello_turn4)
         num_cont_hb(i)=num_conti
       enddo   ! i
+C Loop over all neighbouring boxes
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
 c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
@@ -2838,6 +2958,32 @@ c
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+C Return atom into box, boxxsize is size of box in x dimension
+  164   continue
+        if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
+        if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
+C Condition for being inside the proper box
+        if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
+     &       (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
+        go to 164
+        endif
+  165   continue
+        if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
+        if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
+C Condition for being inside the proper box
+        if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
+     &       (ymedi.lt.((yshift-0.5d0)*boxysize))) then
+        go to 165
+        endif
+  166   continue
+        if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
+        if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+C Condition for being inside the proper box
+        if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
+     &       (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
+        go to 166
+        endif
+
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
@@ -2847,6 +2993,10 @@ c          write (iout,*) i,j,itype(i),itype(j)
         enddo ! j
         num_cont_hb(i)=num_conti
       enddo   ! i
+      enddo   ! zshift
+      enddo   ! yshift
+      enddo   ! xshift
+
 c      write (iout,*) "Number of loop steps in EELEC:",ind
 cd      do i=1,nres
 cd        write (iout,'(i3,3f10.5,5x,3f10.5)') 
@@ -2911,9 +3061,41 @@ c          ind=ind+1
           dx_normj=dc_norm(1,j)
           dy_normj=dc_norm(2,j)
           dz_normj=dc_norm(3,j)
-          xj=c(1,j)+0.5D0*dxj-xmedi
-          yj=c(2,j)+0.5D0*dyj-ymedi
-          zj=c(3,j)+0.5D0*dzj-zmedi
+C          xj=c(1,j)+0.5D0*dxj-xmedi
+C          yj=c(2,j)+0.5D0*dyj-ymedi
+C          zj=c(3,j)+0.5D0*dzj-zmedi
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+C        if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
+  174   continue
+        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+C Condition for being inside the proper box
+        if ((xj.gt.((0.5d0)*boxxsize)).or.
+     &       (xj.lt.((-0.5d0)*boxxsize))) then
+        go to 174
+        endif
+  175   continue
+        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+C Condition for being inside the proper box
+        if ((yj.gt.((0.5d0)*boxysize)).or.
+     &       (yj.lt.((-0.5d0)*boxysize))) then
+        go to 175
+        endif
+  176   continue
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+C Condition for being inside the proper box
+        if ((zj.gt.((0.5d0)*boxzsize)).or.
+     &       (zj.lt.((-0.5d0)*boxzsize))) then
+        go to 176
+        endif
+C        endif !endPBC condintion
+        xj=xj-xmedi
+        yj=yj-ymedi
+        zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
@@ -3278,7 +3460,9 @@ cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &            'eelloc',i,j,eel_loc_ij
-c              write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+           if (eel_loc_ij.ne.0)
+     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
+     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
 
           eel_loc=eel_loc+eel_loc_ij
 C Partial derivatives in virtual-bond dihedral angles gamma
@@ -3667,8 +3851,9 @@ c        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
         call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         eello_turn4=eello_turn4-(s1+s2+s3)
-        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
-     &      'eturn4',i,j,-(s1+s2+s3)
+c             write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
+        if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
+     &      'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
 cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
 cd     &    ' eello_turn4_num',8*eello_turn4_num
 C Derivatives in gamma(i)
@@ -3838,13 +4023,40 @@ C
       r0_scp=4.5d0
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
       do i=iatscp_s,iatscp_e
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-
+C Return atom into box, boxxsize is size of box in x dimension
+  134   continue
+        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
+        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+C Condition for being inside the proper box
+        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
+     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
+        go to 134
+        endif
+  135   continue
+        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
+        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+C Condition for being inside the proper box
+        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
+     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
+        go to 135
+        endif
+  136   continue
+        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
+        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+C Condition for being inside the proper box
+        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
+     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
+        go to 136
+        endif
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
@@ -3855,9 +4067,36 @@ c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
 c         zj=c(3,nres+j)-zi
 C Uncomment following three lines for Ca-p interactions
-          xj=c(1,j)-xi
-          yj=c(2,j)-yi
-          zj=c(3,j)-zi
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+  174   continue
+        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+C Condition for being inside the proper box
+        if ((xj.gt.((0.5d0)*boxxsize)).or.
+     &       (xj.lt.((-0.5d0)*boxxsize))) then
+        go to 174
+        endif
+  175   continue
+        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+C Condition for being inside the proper box
+        if ((yj.gt.((0.5d0)*boxysize)).or.
+     &       (yj.lt.((-0.5d0)*boxysize))) then
+        go to 175
+        endif
+  176   continue
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+C Condition for being inside the proper box
+        if ((zj.gt.((0.5d0)*boxzsize)).or.
+     &       (zj.lt.((-0.5d0)*boxzsize))) then
+        go to 176
+        endif
+          xj=xj-xi
+          yj=yj-yi
+          zj=zj-zi
           rij=xj*xj+yj*yj+zj*zj
           r0ij=r0_scp
           r0ijsq=r0ij*r0ij
@@ -3909,6 +4148,9 @@ cgrad          enddo
 
         enddo ! iint
       enddo ! i
+      enddo !zshift
+      enddo !yshift
+      enddo !xshift
       return
       end
 C-----------------------------------------------------------------------------
@@ -3934,13 +4176,40 @@ C
       evdw2_14=0.0d0
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
       do i=iatscp_s,iatscp_e
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-
+C Return atom into box, boxxsize is size of box in x dimension
+  134   continue
+        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
+        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
+C Condition for being inside the proper box
+        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
+     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
+        go to 134
+        endif
+  135   continue
+        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
+        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
+C Condition for being inside the proper box
+        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
+     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
+        go to 135
+        endif
+  136   continue
+        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
+        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+C Condition for being inside the proper box
+        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
+     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
+        go to 136
+        endif
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
@@ -3951,9 +4220,36 @@ c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
 c         zj=c(3,nres+j)-zi
 C Uncomment following three lines for Ca-p interactions
-          xj=c(1,j)-xi
-          yj=c(2,j)-yi
-          zj=c(3,j)-zi
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+  174   continue
+        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
+        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
+C Condition for being inside the proper box
+        if ((xj.gt.((0.5d0)*boxxsize)).or.
+     &       (xj.lt.((-0.5d0)*boxxsize))) then
+        go to 174
+        endif
+  175   continue
+        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
+        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
+C Condition for being inside the proper box
+        if ((yj.gt.((0.5d0)*boxysize)).or.
+     &       (yj.lt.((-0.5d0)*boxysize))) then
+        go to 175
+        endif
+  176   continue
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+C Condition for being inside the proper box
+        if ((zj.gt.((0.5d0)*boxzsize)).or.
+     &       (zj.lt.((-0.5d0)*boxzsize))) then
+        go to 176
+        endif
+          xj=xj-xi
+          yj=yj-yi
+          zj=zj-zi
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
           fac=rrij**expon2
           e1=fac*fac*aad(itypj,iteli)
@@ -4010,6 +4306,9 @@ cgrad          enddo
 
         enddo ! iint
       enddo ! i
+      enddo !zshift
+      enddo !yshift
+      enddo !xshift
       do i=1,nct
         do j=1,3
           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
@@ -4218,24 +4517,31 @@ c
       estr=0.0d0
       estr1=0.0d0
       do i=ibondp_start,ibondp_end
-        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)
-        else
+        if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
+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)
+c        else
+C       Checking if it involves dummy (NH3+ or COO-) group
+         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+C YES   vbldpDUM is the equlibrium length of spring for Dummy atom
+        diff = vbld(i)-vbldpDUM
+         else
+C NO    vbldp0 is the equlibrium lenght of spring for peptide group
         diff = vbld(i)-vbldp0
-        if (energy_dec) write (iout,*) 
+         endif 
+        if (energy_dec) write (iout,'(a7,i5,4f7.3)') 
      &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
         enddo
 c        write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
-        endif
+c        endif
       enddo
       estr=0.5d0*AKP*estr+estr1
 c
@@ -4373,6 +4679,7 @@ C In following comments this theta will be referred to as t_c.
              bthetk=bthet(k,itype2,ichir21,ichir22)
           endif
          thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
+c         write(iout,*) 'chuj tu', y(k),z(k)
         enddo
         dthett=thet_pred_mean*ssd
         thet_pred_mean=thet_pred_mean*ss+a0thet(it)
@@ -4409,8 +4716,8 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &        E_theta,E_tc)
         endif
         etheta=etheta+ethetai
-        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &      'ebend',i,ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
+     &      'ebend',i,ethetai,theta(i),itype(i)
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
@@ -4431,7 +4738,8 @@ C---------------------------------------------------------------------------
 C Calculate the contributions to both Gaussian lobes.
 C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time)
 C The "polynomial part" of the "standard deviation" of this part of 
-C the distribution.
+C the distributioni.
+        write (iout,*) thetai,thet_pred_mean
         sig=polthet(3,it)
         do j=2,0,-1
           sig=sig*thet_pred_mean+polthet(j,it)
@@ -4461,6 +4769,7 @@ C Following variable is sigma(t_c)**(-2)
         delthe0=thetai-theta0i
         term1=-0.5D0*sigcsq*delthec*delthec
         term2=-0.5D0*sig0inv*delthe0*delthe0
+C        write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0
 C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and
 C NaNs in taking the logarithm. We extract the largest exponent which is added
 C to the energy (this being the log of the distribution) at the end of energy
@@ -4488,6 +4797,7 @@ C Contribution of the bending energy from this theta is just the -log of
 C the sum of the contributions from the two lobes and the pre-exponential
 C factor. Simple enough, isn't it?
         ethetai=(-dlog(termexp)-termm+dlog(termpre))
+C       write (iout,*) 'termexp',termexp,termm,termpre,i
 C NOW the derivatives!!!
 C 6/6/97 Take into account the deformation.
         E_theta=(delthec*sigcsq*term1
@@ -4554,7 +4864,7 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1)) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -5581,8 +5891,12 @@ C Set lprn=.true. for debugging
 c     lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 
-     &       .or. itype(i).eq.ntyp1) cycle
+C ANY TWO ARE DUMMY ATOMS in row CYCLE
+        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
+     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+C For introducing the NH3+ and COO- group please check the etor_d for reference
+C and guidance
         etors_ii=0.0D0
          if (iabs(itype(i)).eq.20) then
          iblock=2
@@ -5683,8 +5997,11 @@ c     lprn=.true.
       etors_d=0.0D0
 c      write(iout,*) "a tu??"
       do i=iphid_start,iphid_end
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
-     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+C ANY TWO ARE DUMMY ATOMS in row CYCLE
+        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or.
+     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))  .or.
+     &      ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -5694,9 +6011,20 @@ c      write(iout,*) "a tu??"
         gloci2=0.0D0
         iblock=1
         if (iabs(itype(i+1)).eq.20) iblock=2
+C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT
+C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO-
+C        if (itype(i+1).eq.ntyp1) iblock=3
+C The problem of NH3+ group can be resolved by adding new parameters please note if there
+C IS or IS NOT need for this
+C IF Yes uncomment below and add to parmread.F appropriate changes and to v1cij and so on
+C        is (itype(i-3).eq.ntyp1) ntblock=2
+C        ntblock is N-terminal blocking group
 
 C Regular cosine and sine terms
         do j=1,ntermd_1(itori,itori1,itori2,iblock)
+C Example of changes for NH3+ blocking group
+C       do j=1,ntermd_1(itori,itori1,itori2,iblock,ntblock)
+C          v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock)
           v1cij=v1c(1,j,itori,itori1,itori2,iblock)
           v1sij=v1s(1,j,itori,itori1,itori2,iblock)
           v2cij=v1c(2,j,itori,itori1,itori2,iblock)

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 6ea0840..bd2165b 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -59,7 +59,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia
 c and Stokes' radii of the peptide group and side chains
 c
 #ifdef CRYST_BOND
-      read (ibond,*) vbldp0,akp,mp,ip,pstok
+      read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
       do i=1,ntyp
         nbondterm(i)=1
         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
@@ -71,7 +71,7 @@ c
         endif
       enddo
 #else
-      read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
+      read (ibond,*) junk,vbldp0,vbldpdum,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)

diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F
index f2daadd..6b385fe 100644 (file)
--- a/source/unres/src_MD-M/readpdb.F
+++ b/source/unres/src_MD-M/readpdb.F
@@ -46,7 +46,8 @@ crc----------------------------------------
           goto 10
         else if (card(:3).eq.'TER') then
 C End current chain
-          ires_old=ires+1 
+          ires_old=ires+2
+          itype(ires_old-1)=ntyp1 
           itype(ires_old)=ntyp1
           ibeg=2
 c          write (iout,*) "Chain ended",ires,ishift,ires_old
@@ -120,13 +121,49 @@ C system
       do i=2,nres-1
 c        write (iout,*) i,itype(i)
         if (itype(i).eq.ntyp1) then
-c          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
-c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
-            dc(j,i)=c(j,i)
-          enddo
-        endif
+         if (itype(i+1).eq.ntyp1) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+           if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+            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
+            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
+             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
+C 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
+            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
 C Calculate the CM of the last side chain.
       if (unres_pdb) then
@@ -150,11 +187,11 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=c(j,nres-2)-c(j,nres-3)
+          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
@@ -181,11 +218,11 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-3.8d0*e2(j)
+            c(j,1)=c(j,2)-1.9d0*e2(j)
           enddo
         else
         do j=1,3
-          dcj=c(j,4)-c(j,3)
+          dcj=(c(j,4)-c(j,3))/2.0
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 265b705..de4b92e 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -213,7 +213,13 @@ cfmc        modecalc=10
       i2ndstr=index(controlcard,'USE_SEC_PRED')
       gradout=index(controlcard,'GRADOUT').gt.0
       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
+C  DISTCHAINMAX become obsolete for periodic boundry condition
       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
+C Reading the dimensions of box in x,y,z coordinates
+      call reada(controlcard,'BOXX',boxxsize,100.0d0)
+      call reada(controlcard,'BOXY',boxysize,100.0d0)
+      call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+
       if (me.eq.king .or. .not.out1file ) 
      &  write (iout,*) "DISTCHAINMAX",distchainmax