eello_turn4=0.0d0
endif
else
-c write (iout,*) "Soft-spheer ELEC potential"
- call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
- & eello_turn4)
+ write (iout,*) "Soft-spheer ELEC potential"
+c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
+c & eello_turn4)
endif
c print *,"Processor",myrank," computed UELEC"
C
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
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)
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
enddo
enddo
#endif
-#undef DEBUG
+c#undef DEBUG
#ifdef DEBUG
write (iout,*) "gloc after reduce"
do i=1,4*nres
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
+C do xshift=-1,1
+C do yshift=-1,1
+C do zshift=-1,1
do i=iatsc_s,iatsc_e
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
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
+C xi=xi+xshift*boxxsize
+C yi=yi+yshift*boxysize
+C zi=zi+zshift*boxzsize
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
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
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+C xj=xj-xi
+C yj=yj-yi
+C 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)
enddo ! j
enddo ! iint
enddo ! i
- enddo ! zshift
- enddo ! yshift
- enddo ! xshift
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
c write (iout,*) "Number of loop steps in EGB:",ind
cccc energy_dec=.false.
return
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)
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)
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)
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)
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)
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)
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 do xshift=-1,1
+C do yshift=-1,1
+C 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
+ & .or. itype(i-1).eq.ntyp1
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
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
+ 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 xmedi=xmedi+xshift*boxxsize
+C ymedi=ymedi+yshift*boxysize
+C zmedi=zmedi+zshift*boxzsize
+
+C Return tom into box, boxxsize is size of box in x dimension
+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 ((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
enddo ! i
- enddo ! zshift
- enddo ! yshift
- enddo ! xshift
+C enddo ! zshift
+C enddo ! yshift
+C enddo ! xshift
c write (iout,*) "Number of loop steps in EELEC:",ind
cd do i=1,nres
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
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ isubchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
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
- zj=zj-zmedi
+C xj=xj-xmedi
+C yj=yj-ymedi
+C zj=zj-zmedi
rij=xj*xj+yj*yj+zj*zj
sss=sscale(sqrt(rij))
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
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
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)
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
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
+C do xshift=-1,1
+C do yshift=-1,1
+C 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)
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 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 ((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
+C xi=xi+xshift*boxxsize
+C yi=yi+yshift*boxysize
+C zi=zi+zshift*boxzsize
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
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
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+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
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+c c endif
+C xj=xj-xi
+C yj=yj-yi
+C zj=zj-zi
rij=xj*xj+yj*yj+zj*zj
r0ij=r0_scp
enddo ! iint
enddo ! i
- enddo !zshift
- enddo !yshift
- enddo !xshift
+C enddo !zshift
+C enddo !yshift
+C enddo !xshift
return
end
C-----------------------------------------------------------------------------
dimension ggg(3)
evdw2=0.0D0
evdw2_14=0.0d0
+c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
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
+C do xshift=-1,1
+C do yshift=-1,1
+C 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))
+ 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 xi=xi+xshift*boxxsize
+c yi=yi+yshift*boxysize
+c zi=zi+zshift*boxzsize
+c print *,xi,yi,zi,'polozenie 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 print *,xi,boxxsize,"pierwszy"
+
+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)
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
+ 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 ((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
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+c if ((zj.gt.((0.5d0)*boxzsize)).or.
+c & (zj.lt.((-0.5d0)*boxzsize))) then
+c go to 176
+c endif
+CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+c print *,xj,yj,zj,'polozenie j'
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+c print *,rrij
sss=sscale(1.0d0/(dsqrt(rrij)))
+c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz'
+c if (sss.eq.0) print *,'czasem jest OK'
+ if (sss.le.0.0d0) cycle
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)
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- endif !endif for sscale cutoff
+c endif !endif for sscale cutoff
enddo ! j
enddo ! iint
enddo ! i
- enddo !zshift
- enddo !yshift
- enddo !xshift
+c enddo !zshift
+c enddo !yshift
+c enddo !xshift
do i=1,nct
do j=1,3
gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)
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)
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
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
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
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
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
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
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...
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
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))
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
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)
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
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))
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),
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)
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,
projects / unres.git / blobdiff