X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=346a3c3200b1acae0d76e979343576bac8e5bb39;hb=17450e1f47ecb1b7f7cb5c9008ebc14164ca69aa;hp=8f89ae98b154ad61c4d720fe2ced06fbc053a9e5;hpb=d609f4b92f230a11e364e22af6514fe261b5c89d;p=unres.git 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 8f89ae9..346a3c3 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -711,7 +711,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -720,7 +720,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -740,7 +740,7 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -749,7 +749,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -1436,30 +1436,36 @@ C we have the original box) 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 134 continue +c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize +c 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 if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 134 +c endif +c 135 continue +c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize +c 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-boxzsize - if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. +c & (yi.lt.((yshift-0.5d0)*boxysize))) then +c go to 135 +c endif +c 136 continue +c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize +c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize 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 +c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 136 +c endif + 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) @@ -1505,31 +1511,36 @@ c alf12=0.0D0 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 137 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c 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 if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 137 +c endif +c 138 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c 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-boxzsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 138 +c endif +c 139 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize 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 - +c if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 139 +c endif + 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 dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -2404,13 +2415,13 @@ c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itortyp(itype(i-2)) else - iti=ntortyp+1 + iti=ntortyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then iti1 = itortyp(itype(i-1)) else - iti1=ntortyp+1 + iti1=ntortyp endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) @@ -2451,10 +2462,10 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (itype(i-1).le.ntyp) then iti1 = itortyp(itype(i-1)) else - iti1=ntortyp+1 + iti1=ntortyp endif else - iti1=ntortyp+1 + iti1=ntortyp endif do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) @@ -2869,7 +2880,11 @@ 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 + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2880,30 +2895,36 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition 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-boxzsize - if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +c 184 continue +c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize +c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+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 +c if ((xmedi.gt.((0.5d0)*boxxsize)).or. +c & (xmedi.lt.((-0.5d0)*boxxsize))) then +c go to 184 +c endif +c 185 continue +c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize +c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize +cC Condition for being inside the proper box +c if ((ymedi.gt.((0.5d0)*boxysize)).or. +c & (ymedi.lt.((-0.5d0)*boxysize))) then +c go to 185 +c endif +c 186 continue +c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize +c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +cC Condition for being inside the proper box +c if ((zmedi.gt.((0.5d0)*boxzsize)).or. +c & (zmedi.lt.((-0.5d0)*boxzsize))) then +c go to 186 +c endif + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2912,7 +2933,11 @@ C Condition for being inside the proper box do i=iturn4_start,iturn4_end if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 - & .or. itype(i+4).eq.ntyp1) cycle + & .or. itype(i+4).eq.ntyp1 + & .or. itype(i+5).eq.ntyp1 + & .or. itype(i).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2923,30 +2948,36 @@ C Condition for being inside the proper box 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 194 continue +c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize +c 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 if ((xmedi.gt.((0.5d0)*boxxsize)).or. +c & (xmedi.lt.((-0.5d0)*boxxsize))) then +c go to 194 +c endif +c 195 continue +c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize +c 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-boxzsize - if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +c if ((ymedi.gt.((0.5d0)*boxysize)).or. +c & (ymedi.lt.((-0.5d0)*boxysize))) then +c go to 195 +c endif +c 196 continue +c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize +c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize 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 +c if ((zmedi.gt.((0.5d0)*boxzsize)).or. +c & (zmedi.lt.((-0.5d0)*boxzsize))) then +c go to 196 +c endif + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=num_cont_hb(i) call eelecij(i,i+3,ees,evdw1,eel_loc) @@ -2962,7 +2993,10 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i-1).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2972,37 +3006,47 @@ c xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize + 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 164 continue +c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize +c 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-boxzsize - if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 164 +c endif +c 165 continue +c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize +c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize 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 if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or. +c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then +c go to 165 +c endif +c 166 continue +c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize +c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +cC Condition for being inside the proper box +c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 166 +c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) c write (iout,*) i,j,itype(i),itype(j) - if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle + if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + & .or.itype(j-1).eq.ntyp1 + &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -3082,31 +3126,38 @@ 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 + 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 + 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 174 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c 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 if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 174 +c endif +c 175 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c 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-boxzsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 175 +c endif +c 176 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize 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 if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 176 +c endif C endif !endPBC condintion xj=xj-xmedi yj=yj-ymedi @@ -3186,9 +3237,15 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + if (sss.gt.0.0) then + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + else + ggg(1)=0.0 + ggg(2)=0.0 + ggg(3)=0.0 + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -3212,7 +3269,7 @@ C MARYSIA facvdw=(ev1+evdwij)*sss facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel)+sssgrad*rmij*evdwij + fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -3241,9 +3298,9 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end - ggg(1)=facvdw*xj*sss - ggg(2)=facvdw*yj*sss - ggg(3)=facvdw*zj*sss + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -3479,7 +3536,8 @@ cgrad endif C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) -cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +c write (iout,*) 'i',i,' j',j,itype(i),itype(j), +c & ' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij @@ -4056,30 +4114,37 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e 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 134 continue +c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize +c 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-boxzsize - if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 134 +c endif +c 135 continue +c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize +c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize 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 +c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. +c & (yi.lt.((yshift-0.5d0)*boxysize))) then +c go to 135 +c c endif +c 136 continue +c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize +c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +cC Condition for being inside the proper box +c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 136 +c endif + 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 + do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -4093,30 +4158,36 @@ C Uncomment following three lines for Ca-p interactions 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 174 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize 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-boxzsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +c if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 174 +c endif +c 175 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +cC Condition for being inside the proper box +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 175 +c endif +c 176 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize 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 if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 176 + 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 +c c endif xj=xj-xi yj=yj-yi zj=zj-zi @@ -4210,31 +4281,38 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e 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)) + 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 + 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 134 continue +c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize +c 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 if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. +c & (xi.lt.((xshift-0.5d0)*boxxsize))) then +c go to 134 +c endif +c 135 continue +c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize +c 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-boxzsize - if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize +c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. +c & (yi.lt.((yshift-0.5d0)*boxysize))) then +c go to 135 +c endif +c 136 continue +c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize +c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize 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 +c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. +c & (zi.lt.((zshift-0.5d0)*boxzsize))) then +c go to 136 +c endif do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -4248,30 +4326,36 @@ C Uncomment following three lines for Ca-p interactions 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 + 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 +c 174 continue +c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize +c 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-boxzsize - if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize +c if ((xj.gt.((0.5d0)*boxxsize)).or. +c & (xj.lt.((-0.5d0)*boxxsize))) then +c go to 174 +c endif +c 175 continue +c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize +c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +cC Condition for being inside the proper box +c if ((yj.gt.((0.5d0)*boxysize)).or. +c & (yj.lt.((-0.5d0)*boxysize))) then +c go to 175 +c endif +c 176 continue +c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize +c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize 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 if ((zj.gt.((0.5d0)*boxzsize)).or. +c & (zj.lt.((-0.5d0)*boxzsize))) then +c go to 176 +c endif xj=xj-xi yj=yj-yi zj=zj-zi @@ -4652,9 +4736,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end - print *,i,itype(i-1),itype(i),itype(i-2) - if (itype(i-1).eq.ntyp1) cycle - print *,'wchodze',itype(i-1) + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) @@ -4671,7 +4754,7 @@ C Zero the energy function and its derivative at 0 or pi. ichir22=isign(1,itype(i)) endif - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4684,7 +4767,7 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4692,8 +4775,8 @@ C Zero the energy function and its derivative at 0 or pi. z(1)=cos(phii1) #else phii1=phi(i+1) - z(1)=dcos(phii1) #endif + z(1)=dcos(phii1) z(2)=dsin(phii1) else z(1)=0.0D0 @@ -4897,7 +4980,10 @@ C etheta=0.0D0 do i=ithet_start,ithet_end c print *,i,itype(i-1),itype(i),itype(i-2) - if ((itype(i-1).eq.ntyp1)) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF + if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -4909,7 +4995,7 @@ c print *,i,itype(i-1),itype(i),itype(i-2) coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4930,7 +5016,7 @@ C propagation of chirality for glycine type sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -5828,9 +5914,9 @@ c lprn=.true. do i=iphi_start,iphi_end etors_ii=0.0D0 if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Proline-Proline pair is a special case... @@ -5928,8 +6014,9 @@ C ANY TWO ARE DUMMY ATOMS in row CYCLE c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle - if ((itype(i-3).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. - & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle + if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF C For introducing the NH3+ and COO- group please check the etor_d for reference C and guidance etors_ii=0.0D0 @@ -6033,10 +6120,14 @@ c lprn=.true. c write(iout,*) "a tu??" do i=iphid_start,iphid_end 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 +C if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +C & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or. +C & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1)) .or. +C & ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle + if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. + & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or. + & (itype(i+1).eq.ntyp1)) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -6046,6 +6137,7 @@ C ANY TWO ARE DUMMY ATOMS in row CYCLE gloci2=0.0D0 iblock=1 if (iabs(itype(i+1)).eq.20) iblock=2 +C Iblock=2 Proline type 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 @@ -7174,7 +7266,7 @@ C--------------------------------------------------------------------------- if (j.lt.nres-1) then itj1 = itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -7264,14 +7356,14 @@ C parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then iti=itortyp(itype(i)) else - iti=ntortyp+1 + iti=ntortyp endif itk1=itortyp(itype(k+1)) itj=itortyp(itype(j)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -7417,7 +7509,7 @@ C Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then iti=itortyp(itype(i)) else - iti=ntortyp+1 + iti=ntortyp endif itk1=itortyp(itype(k+1)) itl=itortyp(itype(l)) @@ -7425,7 +7517,7 @@ C Antiparallel orientation of the two CA-CA-CA frames. if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif C A2 kernel(j-1)T A1T call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), @@ -8579,14 +8671,14 @@ C energy moment and not to the cluster cumulant. if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif itk=itortyp(itype(k)) itk1=itortyp(itype(k+1)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) @@ -8698,19 +8790,19 @@ cd write (2,*) 'eello_graph4: wturn6',wturn6 if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif itk=itortyp(itype(k)) if (k.lt.nres-1) then itk1=itortyp(itype(k+1)) else - itk1=ntortyp+1 + itk1=ntortyp endif itl=itortyp(itype(l)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,