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=e56e104ce46f8f20ca3675d59eb08e552b98131f;hb=892c42fe1b098aac20f41e47faf97b1b2479b647;hp=fa388be8119f2f3171e65d7b5698e44d8f86010c;hpb=24267f5c78a1b1b5b8a5eaec6ec89a085630502e;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 fa388be..e56e104 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -24,6 +24,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.MD' include 'COMMON.CONTROL' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -296,6 +297,8 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor +c Here are the energies showed per procesor if the are more processors +c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) c print *," Processor",myrank," left SUM_ENERGY" @@ -387,14 +390,14 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor @@ -1410,7 +1413,9 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' logical lprn + integer xshift,yshift,zshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -1418,6 +1423,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 @@ -1425,6 +1435,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-boxzsize + 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 + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1465,17 +1501,52 @@ 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-boxzsize + 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 + 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) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) + sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) + +c write (iout,'(a7,4f8.3)') +c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb + if (sss.gt.0.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -1504,7 +1575,7 @@ c--------------------------------------------------------------- c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij + evdw=evdw+evdwij*sss if (lprn) then sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -1524,6 +1595,9 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac +c print '(2i4,6f8.4)',i,j,sss,sssgrad* +c & evdwij,fac,sigma(itypi,itypj),expon + fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij c fac=0.0d0 C Calculate the radial part of the gradient gg(1)=xj*fac @@ -1531,9 +1605,13 @@ C Calculate the radial part of the gradient gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad + endif 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 @@ -1741,6 +1819,7 @@ C---------------------------------------------------------------------------- include 'COMMON.CALC' include 'COMMON.IOUNITS' double precision dcosom1(3),dcosom2(3) +cc print *,'sss=',sss eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 @@ -1759,16 +1838,16 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 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 enddo c 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 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 c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) @@ -2325,13 +2404,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) @@ -2372,10 +2451,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) @@ -2704,6 +2783,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -2786,9 +2866,14 @@ 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 + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 +c & .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) @@ -2798,6 +2883,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-boxzsize + 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 186 + endif num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2806,7 +2916,9 @@ C 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 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2816,17 +2928,49 @@ 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-boxzsize + 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 + 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 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 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2836,15 +2980,47 @@ 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-boxzsize + if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize +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) 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 + &) cycle call eelecij(i,j,ees,evdw1,eel_loc) 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)') @@ -2875,6 +3051,7 @@ C------------------------------------------------------------------------------- include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -2909,10 +3086,46 @@ 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-boxzsize + 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 endif !endPBC condintion + xj=xj-xmedi + yj=yj-ymedi + zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj + + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) +c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij @@ -2928,14 +3141,15 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions ev2=bbb*r6ij fac3=ael6i*r6ij fac4=ael3i*r3ij - evdwij=ev1+ev2 + evdwij=(ev1+ev2) el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac - eesij=el1+el2 +C MARYSIA + eesij=(el1+el2) C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) ees=ees+eesij - evdw1=evdw1+evdwij + evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, @@ -2952,7 +3166,7 @@ C C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij) + facvdw=-6*rrmij*(ev1+evdwij)*sss facel=-3*rrmij*(el1+eesij) fac1=fac erij(1)=xj*rmij @@ -2982,9 +3196,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 @@ -3004,8 +3224,9 @@ cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else - facvdw=ev1+evdwij - facel=el1+eesij +C MARYSIA + facvdw=(ev1+evdwij)*sss + facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij @@ -3036,9 +3257,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 - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + 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) @@ -3077,14 +3298,16 @@ cgrad enddo cgrad enddo do k=1,3 gelc(k,i)=gelc(k,i) - & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) gelc(k,j)=gelc(k,j) - & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C MARYSIA +c endif !sscale IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN @@ -3272,11 +3495,14 @@ 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 -c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) +c if (eel_loc_ij.ne.0) +c & write (iout,'(a4,2i4,8f9.5)')'chuj', +c & 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 @@ -3304,14 +3530,14 @@ cgrad enddo cgrad enddo C Remaining derivatives of eello do l=1,3 - gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ - & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4) - gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+ - & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4) - gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+ - & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4) - gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+ - & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4) + gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ + & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) + gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ + & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) + gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ + & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) + gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ + & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) enddo ENDIF C Change 12/26/95 to calculate four-body contributions to H-bonding energy @@ -3665,8 +3891,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) @@ -3836,13 +4063,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-boxzsize + 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 do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -3853,10 +4107,38 @@ 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-boxzsize + 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 + xj=xj-xi + yj=yj-yi + zj=zj-zi rij=xj*xj+yj*yj+zj*zj + r0ij=r0_scp r0ijsq=r0ij*r0ij if (rij.lt.r0ijsq) then @@ -3907,6 +4189,9 @@ cgrad enddo enddo ! iint enddo ! i + enddo !zshift + enddo !yshift + enddo !xshift return end C----------------------------------------------------------------------------- @@ -3927,18 +4212,46 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' dimension ggg(3) evdw2=0.0D0 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-boxzsize + 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 do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -3949,27 +4262,58 @@ 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-boxzsize + 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 + xj=xj-xi + yj=yj-yi + zj=zj-zi rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + sss=sscale(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) + if (sss.gt.0.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+e1+e2 + evdw2_14=evdw2_14+(e1+e2)*sss endif evdwij=e1+e2 - evdw2=evdw2+evdwij + evdw2=evdw2+evdwij*sss if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C - fac=-(evdwij+e1)*rrij + fac=-(evdwij+e1)*rrij*sss + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -4004,10 +4348,14 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - enddo + endif !endif for sscale cutoff + enddo ! j 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) @@ -4163,7 +4511,7 @@ c dscj_inv=dsc_inv(itypj) cosphi=om12-om1*om2 eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) & +akct*deltad*deltat12 - & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi + & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, c & " deltat12",deltat12," eij",eij @@ -4216,24 +4564,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 @@ -4314,7 +4669,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end - 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 Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) @@ -4331,7 +4687,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 @@ -4344,7 +4700,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 @@ -4352,8 +4708,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 @@ -4371,6 +4727,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) @@ -4407,11 +4764,11 @@ 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) + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo C Ufff.... We've done all this!!! return @@ -4429,7 +4786,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. +ccc write (iout,*) thetai,thet_pred_mean sig=polthet(3,it) do j=2,0,-1 sig=sig*thet_pred_mean+polthet(j,it) @@ -4459,6 +4817,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 @@ -4486,6 +4845,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 @@ -4552,7 +4912,11 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle +c print *,i,itype(i-1),itype(i),itype(i-2) + 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 @@ -4564,7 +4928,7 @@ C 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 @@ -4585,7 +4949,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 @@ -4717,7 +5081,7 @@ c lprn1=.false. etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=wang*dethetai+ gloc(nphi+i-2,icg) enddo return end @@ -5114,7 +5478,7 @@ c write (2,*) "xx",xx," yy",yy," zz",zz Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@ -5483,9 +5847,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... @@ -5579,8 +5943,15 @@ 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 +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-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 if (iabs(itype(i)).eq.20) then iblock=2 @@ -5681,8 +6052,15 @@ 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 +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)) @@ -5692,9 +6070,21 @@ c write(iout,*) "a tu??" 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 +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) @@ -6809,7 +7199,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) @@ -6899,14 +7289,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 @@ -7052,7 +7442,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)) @@ -7060,7 +7450,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), @@ -8214,14 +8604,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) @@ -8333,19 +8723,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,