source/unres/src-5hdiag-tmp/edif

   1 27d26
   2 <       include 'COMMON.SPLITELE'
   3 102d100
   4 < C      print *,ipot
   5 116d113
   6 < C      print *,"bylem w egb"
   7 127,131d123
   8 < cmc
   9 < cmc Sep-06: egb takes care of dynamic ss bonds too
  10 < cmc
  11 < c      if (dyn_ss) call dyn_set_nss
  12 <
  13 140,149d131
  14 < C Introduction of shielding effect first for each peptide group
  15 < C the shielding factor is set this factor is describing how each
  16 < C peptide group is shielded by side-chains
  17 < C the matrix - shield_fac(i) the i index describe the ith between i and i+1
  18 < C      write (iout,*) "shield_mode",shield_mode
  19 <       if (shield_mode.eq.1) then
  20 <        call set_shield_fac
  21 <       else if  (shield_mode.eq.2) then
  22 <        call set_shield_fac2
  23 <       endif
  24 172,174c154,156
  25 <         write (iout,*) "Soft-spheer ELEC potential"
  26 < c        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
  27 < c     &   eello_turn4)
  28 ---
  29 > c        write (iout,*) "Soft-spheer ELEC potential"
  30 >         call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
  31 >      &   eello_turn4)
  32 206,213c188
  33 <        if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
  34 <         call ebend(ebe,ethetacnstr)
  35 <         endif
  36 < C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
  37 < C energy function
  38 <        if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
  39 <          call ebend_kcc(ebe,ethetacnstr)
  40 <         endif
  41 ---
  42 >         call ebend(ebe)
  43 216d190
  44 <         ethetacnstr=0
  45 222d195
  46 < C      print *,"TU DOCHODZE?"
  47 229d201
  48 < C      print *,"tor",tor_mode
  49 231d202
  50 <        if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
  51 233,238d203
  52 <        endif
  53 < C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
  54 < C energy function
  55 <        if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
  56 <        call etor_kcc(etors,edihcnstr)
  57 <        endif
  58 247c212
  59 <       if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then
  60 ---
  61 >       if (wtor_d.gt.0) then
  62 261d225
  63 < C      print *,"PRZED MULIt"
  64 294,306d257
  65 < C 01/27/2015 added by adasko
  66 < C the energy component below is energy transfer into lipid environment
  67 < C based on partition function
  68 < C      print *,"przed lipidami"
  69 <       if (wliptran.gt.0) then
  70 <         call Eliptransfer(eliptran)
  71 <       endif
  72 < C      print *,"za lipidami"
  73 <       if (AFMlog.gt.0) then
  74 <         call AFMforce(Eafmforce)
  75 <       else if (selfguide.gt.0) then
  76 <         call AFMvel(Eafmforce)
  77 <       endif
  78 348,350d298
  79 <       energia(22)=eliptran
  80 <       energia(23)=Eafmforce
  81 <       energia(24)=ethetacnstr
  82 355d302
  83 <       if (dyn_ss) call dyn_set_nss
  84 442,444d388
  85 <       eliptran=energia(22)
  86 <       Eafmforce=energia(23)
  87 <       ethetacnstr=energia(24)
  88 451,452c395
  89 <      & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
  90 <      & +ethetacnstr
  91 ---
  92 >      & +wbond*estr+Uconst+wsccor*esccor
  93 459,461c402
  94 <      & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
  95 <      & +Eafmforce
  96 <      & +ethetacnstr
  97 ---
  98 >      & +wbond*estr+Uconst+wsccor*esccor
  99 496a438,439
 100 >       double precision gradbufc(3,maxres),gradbufx(3,maxres),
 101 >      &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
 102 498,500d440
 103 <       double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres),
 104 <      & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres)
 105 <      & ,gloc_scbuf(3,-1:maxres)
 106 553c493
 107 <       do i=0,nct
 108 ---
 109 >       do i=1,nct
 110 564,572d503
 111 <      &                +wliptran*gliptranc(j,i)
 112 <      &                +gradafm(j,i)
 113 <      &                 +welec*gshieldc(j,i)
 114 <      &                 +wcorr*gshieldc_ec(j,i)
 115 <      &                 +wturn3*gshieldc_t3(j,i)
 116 <      &                 +wturn4*gshieldc_t4(j,i)
 117 <      &                 +wel_loc*gshieldc_ll(j,i)
 118 <
 119 <
 120 576c507
 121 <       do i=0,nct
 122 ---
 123 >       do i=1,nct
 124 588,595d518
 125 <      &                +wliptran*gliptranc(j,i)
 126 <      &                +gradafm(j,i)
 127 <      &                 +welec*gshieldc(j,i)
 128 <      &                 +wcorr*gshieldc_ec(j,i)
 129 <      &                 +wturn4*gshieldc_t4(j,i)
 130 <      &                 +wel_loc*gshieldc_ll(j,i)
 131 <
 132 <
 133 609c532
 134 <       do i=0,nres
 135 ---
 136 >       do i=1,nres
 137 652c575
 138 <       do i=nres-2,-1,-1
 139 ---
 140 >       do i=nres-2,nnt,-1
 141 673c596
 142 <       do i=-1,nres
 143 ---
 144 >       do i=1,nres
 145 682c605
 146 <       do i=nres-2,-1,-1
 147 ---
 148 >       do i=nres-2,nnt,-1
 149 710c633
 150 <       do i=-1,nct
 151 ---
 152 >       do i=1,nct
 153 713,719d635
 154 < C          print *,gradbufc(1,13)
 155 < C          print *,welec*gelc(1,13)
 156 < C          print *,wel_loc*gel_loc(1,13)
 157 < C          print *,0.5d0*(wscp*gvdwc_scpp(1,13))
 158 < C          print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13)
 159 < C          print *,wel_loc*gel_loc_long(1,13)
 160 < C          print *,gradafm(1,13),"AFM"
 161 738,755d653
 162 <      &               +wliptran*gliptranc(j,i)
 163 <      &                +gradafm(j,i)
 164 <      &                 +welec*gshieldc(j,i)
 165 <      &                 +welec*gshieldc_loc(j,i)
 166 <      &                 +wcorr*gshieldc_ec(j,i)
 167 <      &                 +wcorr*gshieldc_loc_ec(j,i)
 168 <      &                 +wturn3*gshieldc_t3(j,i)
 169 <      &                 +wturn3*gshieldc_loc_t3(j,i)
 170 <      &                 +wturn4*gshieldc_t4(j,i)
 171 <      &                 +wturn4*gshieldc_loc_t4(j,i)
 172 <      &                 +wel_loc*gshieldc_ll(j,i)
 173 <      &                 +wel_loc*gshieldc_loc_ll(j,i)
 174 <
 175 <
 176 <
 177 <
 178 <
 179 <
 180 760c658
 181 <      &                welec*gelc_long(j,i)+
 182 ---
 183 >      &                welec*gelc_long(j,i)
 184 775,791d672
 185 <      &               +wliptran*gliptranc(j,i)
 186 <      &                +gradafm(j,i)
 187 <      &                 +welec*gshieldc(j,i)
 188 <      &                 +welec*gshieldc_loc(j,i)
 189 <      &                 +wcorr*gshieldc_ec(j,i)
 190 <      &                 +wcorr*gshieldc_loc_ec(j,i)
 191 <      &                 +wturn3*gshieldc_t3(j,i)
 192 <      &                 +wturn3*gshieldc_loc_t3(j,i)
 193 <      &                 +wturn4*gshieldc_t4(j,i)
 194 <      &                 +wturn4*gshieldc_loc_t4(j,i)
 195 <      &                 +wel_loc*gshieldc_ll(j,i)
 196 <      &                 +wel_loc*gshieldc_loc_ll(j,i)
 197 <
 198 <
 199 <
 200 <
 201 <
 202 798,806d678
 203 <      &                 +wliptran*gliptranx(j,i)
 204 <      &                 +welec*gshieldx(j,i)
 205 <      &                 +wcorr*gshieldx_ec(j,i)
 206 <      &                 +wturn3*gshieldx_t3(j,i)
 207 <      &                 +wturn4*gshieldx_t4(j,i)
 208 <      &                 +wel_loc*gshieldx_ll(j,i)
 209 <
 210 <
 211 <
 212 841d712
 213 < c#define DEBUG
 214 850d720
 215 < c#undef DEBUG
 216 870d739
 217 < c#define DEBUG
 218 879d747
 219 < c#undef DEBUG
 220 1009d876
 221 <       include 'COMMON.CONTROL'
 222 1045,1049d911
 223 <       if (shield_mode.gt.0) then
 224 <        wscp=weights(2)*fact
 225 <        wsc=weights(1)*fact
 226 <        wvdwpp=weights(16)*fact
 227 <       endif
 228 1101,1103d962
 229 <       eliptran=energia(22)
 230 <       Eafmforce=energia(23)
 231 <       ethetacnstr=energia(24)
 232 1110,1112c969,971
 233 <      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
 234 <      &  ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
 235 <      &  etot
 236 ---
 237 >      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
 238 >      &  edihcnstr,ebr*nss,
 239 >      &  Uconst,etot
 240 1134d992
 241 <      & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
 242 1137,1138d994
 243 <      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
 244 <      & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
 245 1140d995
 246 <
 247 1148,1149c1003
 248 <      &  ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
 249 <      &  etot
 250 ---
 251 >      &  ebr*nss,Uconst,etot
 252 1170d1023
 253 <      & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
 254 1173,1174d1025
 255 <      & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
 256 <      & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
 257 1229,1231c1080,1081
 258 < C have you changed here?
 259 <             e1=fac*fac*aa
 260 <             e2=fac*bb
 261 ---
 262 >             e1=fac*fac*aa(itypi,itypj)
 263 >             e2=fac*bb(itypi,itypj)
 264 1236c1086
 265 < cd   &        restyp(itypi),i,restyp(itypj),j,a(itypi,itypj),
 266 ---
 267 > cd   &        restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 268 1380,1382c1230,1231
 269 < C have you changed here?
 270 <             e1=fac*fac*aa
 271 <             e2=fac*bb
 272 ---
 273 >             e1=fac*fac*aa(itypi,itypj)
 274 >             e2=fac*bb(itypi,itypj)
 275 1508d1356
 276 < C have you changed here?
 277 1510,1511c1358,1359
 278 <             e1=fac*fac*aa
 279 <             e2=fac*bb
 280 ---
 281 >             e1=fac*fac*aa(itypi,itypj)
 282 >             e2=fac*bb(itypi,itypj)
 283 1518,1519c1366,1367
 284 <             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 285 <             epsi=bb**2/aa
 286 ---
 287 >             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 288 >             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 289 1563,1564d1410
 290 <       include 'COMMON.SPLITELE'
 291 <       include 'COMMON.SBRIDGE'
 292 1566,1567d1411
 293 <       integer xshift,yshift,zshift
 294 <
 295 1570c1414
 296 < C      print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 297 ---
 298 > c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 299 1575,1579d1418
 300 < C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0
 301 < C we have the original box)
 302 < C      do xshift=-1,1
 303 < C      do yshift=-1,1
 304 < C      do zshift=-1,1
 305 1587,1650d1425
 306 < C Return atom into box, boxxsize is size of box in x dimension
 307 < c  134   continue
 308 < c        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
 309 < c        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
 310 < C Condition for being inside the proper box
 311 < c        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
 312 < c     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
 313 < c        go to 134
 314 < c        endif
 315 < c  135   continue
 316 < c        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
 317 < c        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
 318 < C Condition for being inside the proper box
 319 < c        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
 320 < c     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
 321 < c        go to 135
 322 < c        endif
 323 < c  136   continue
 324 < c        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
 325 < c        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
 326 < C Condition for being inside the proper box
 327 < c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
 328 < c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
 329 < c        go to 136
 330 < c        endif
 331 <           xi=mod(xi,boxxsize)
 332 <           if (xi.lt.0) xi=xi+boxxsize
 333 <           yi=mod(yi,boxysize)
 334 <           if (yi.lt.0) yi=yi+boxysize
 335 <           zi=mod(zi,boxzsize)
 336 <           if (zi.lt.0) zi=zi+boxzsize
 337 < C define scaling factor for lipids
 338 <
 339 < C        if (positi.le.0) positi=positi+boxzsize
 340 < C        print *,i
 341 < C first for peptide groups
 342 < c for each residue check if it is in lipid or lipid water border area
 343 <        if ((zi.gt.bordlipbot)
 344 <      &.and.(zi.lt.bordliptop)) then
 345 < C the energy transfer exist
 346 <         if (zi.lt.buflipbot) then
 347 < C what fraction I am in
 348 <          fracinbuf=1.0d0-
 349 <      &        ((zi-bordlipbot)/lipbufthick)
 350 < C lipbufthick is thickenes of lipid buffore
 351 <          sslipi=sscalelip(fracinbuf)
 352 <          ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
 353 <         elseif (zi.gt.bufliptop) then
 354 <          fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
 355 <          sslipi=sscalelip(fracinbuf)
 356 <          ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
 357 <         else
 358 <          sslipi=1.0d0
 359 <          ssgradlipi=0.0
 360 <         endif
 361 <        else
 362 <          sslipi=0.0d0
 363 <          ssgradlipi=0.0
 364 <        endif
 365 <
 366 < C          xi=xi+xshift*boxxsize
 367 < C          yi=yi+yshift*boxysize
 368 < C          zi=zi+zshift*boxzsize
 369 <
 370 1663,1694d1437
 371 <             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
 372 <
 373 < c              write(iout,*) "PRZED ZWYKLE", evdwij
 374 <               call dyn_ssbond_ene(i,j,evdwij)
 375 < c              write(iout,*) "PO ZWYKLE", evdwij
 376 <
 377 <               evdw=evdw+evdwij
 378 <               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
 379 <      &                        'evdw',i,j,evdwij,' ss'
 380 < C triple bond artifac removal
 381 <              do k=j+1,iend(i,iint)
 382 < C search over all next residues
 383 <               if (dyn_ss_mask(k)) then
 384 < C check if they are cysteins
 385 < C              write(iout,*) 'k=',k
 386 <
 387 < c              write(iout,*) "PRZED TRI", evdwij
 388 <                evdwij_przed_tri=evdwij
 389 <               call triple_ssbond_ene(i,j,k,evdwij)
 390 < c               if(evdwij_przed_tri.ne.evdwij) then
 391 < c                 write (iout,*) "TRI:", evdwij, evdwij_przed_tri
 392 < c               endif
 393 <
 394 < c              write(iout,*) "PO TRI", evdwij
 395 < C call the energy function that removes the artifical triple disulfide
 396 < C bond the soubroutine is located in ssMD.F
 397 <               evdw=evdw+evdwij
 398 <               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
 399 <      &                        'evdw',i,j,evdwij,'tss'
 400 <               endif!dyn_ss_mask(k)
 401 <              enddo! k
 402 <             ELSE
 403 1723,1817c1466,1468
 404 <             xj=c(1,nres+j)
 405 <             yj=c(2,nres+j)
 406 <             zj=c(3,nres+j)
 407 < C Return atom J into box the original box
 408 < c  137   continue
 409 < c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
 410 < c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
 411 < C Condition for being inside the proper box
 412 < c        if ((xj.gt.((0.5d0)*boxxsize)).or.
 413 < c     &       (xj.lt.((-0.5d0)*boxxsize))) then
 414 < c        go to 137
 415 < c        endif
 416 < c  138   continue
 417 < c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
 418 < c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
 419 < C Condition for being inside the proper box
 420 < c        if ((yj.gt.((0.5d0)*boxysize)).or.
 421 < c     &       (yj.lt.((-0.5d0)*boxysize))) then
 422 < c        go to 138
 423 < c        endif
 424 < c  139   continue
 425 < c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
 426 < c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
 427 < C Condition for being inside the proper box
 428 < c        if ((zj.gt.((0.5d0)*boxzsize)).or.
 429 < c     &       (zj.lt.((-0.5d0)*boxzsize))) then
 430 < c        go to 139
 431 < c        endif
 432 <           xj=mod(xj,boxxsize)
 433 <           if (xj.lt.0) xj=xj+boxxsize
 434 <           yj=mod(yj,boxysize)
 435 <           if (yj.lt.0) yj=yj+boxysize
 436 <           zj=mod(zj,boxzsize)
 437 <           if (zj.lt.0) zj=zj+boxzsize
 438 <        if ((zj.gt.bordlipbot)
 439 <      &.and.(zj.lt.bordliptop)) then
 440 < C the energy transfer exist
 441 <         if (zj.lt.buflipbot) then
 442 < C what fraction I am in
 443 <          fracinbuf=1.0d0-
 444 <      &        ((zj-bordlipbot)/lipbufthick)
 445 < C lipbufthick is thickenes of lipid buffore
 446 <          sslipj=sscalelip(fracinbuf)
 447 <          ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
 448 <         elseif (zj.gt.bufliptop) then
 449 <          fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
 450 <          sslipj=sscalelip(fracinbuf)
 451 <          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
 452 <         else
 453 <          sslipj=1.0d0
 454 <          ssgradlipj=0.0
 455 <         endif
 456 <        else
 457 <          sslipj=0.0d0
 458 <          ssgradlipj=0.0
 459 <        endif
 460 <       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 461 <      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 462 <       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 463 <      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 464 < C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
 465 < C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 466 < C      if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
 467 < C      print *,sslipi,sslipj,bordlipbot,zi,zj
 468 <       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 469 <       xj_safe=xj
 470 <       yj_safe=yj
 471 <       zj_safe=zj
 472 <       subchap=0
 473 <       do xshift=-1,1
 474 <       do yshift=-1,1
 475 <       do zshift=-1,1
 476 <           xj=xj_safe+xshift*boxxsize
 477 <           yj=yj_safe+yshift*boxysize
 478 <           zj=zj_safe+zshift*boxzsize
 479 <           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 480 <           if(dist_temp.lt.dist_init) then
 481 <             dist_init=dist_temp
 482 <             xj_temp=xj
 483 <             yj_temp=yj
 484 <             zj_temp=zj
 485 <             subchap=1
 486 <           endif
 487 <        enddo
 488 <        enddo
 489 <        enddo
 490 <        if (subchap.eq.1) then
 491 <           xj=xj_temp-xi
 492 <           yj=yj_temp-yi
 493 <           zj=zj_temp-zi
 494 <        else
 495 <           xj=xj_safe-xi
 496 <           yj=yj_safe-yi
 497 <           zj=zj_safe-zi
 498 <        endif
 499 ---
 500 >             xj=c(1,nres+j)-xi
 501 >             yj=c(2,nres+j)-yi
 502 >             zj=c(3,nres+j)-zi
 503 1821,1823d1471
 504 < C            xj=xj-xi
 505 < C            yj=yj-yi
 506 < C            zj=zj-zi
 507 1829,1834d1476
 508 <             sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
 509 <             sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
 510 <
 511 < c            write (iout,'(a7,4f8.3)')
 512 < c    &      "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
 513 <             if (sss.gt.0.0d0) then
 514 1855,1859c1497,1498
 515 < C here to start with
 516 < C            if (c(i,3).gt.
 517 <             faclip=fac
 518 <             e1=fac*fac*aa
 519 <             e2=fac*bb
 520 ---
 521 >             e1=fac*fac*aa(itypi,itypj)
 522 >             e2=fac*bb(itypi,itypj)
 523 1863,1865d1501
 524 < C       write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
 525 < C     &((sslipi+sslipj)/2.0d0+
 526 < C     &(2.0d0-sslipi-sslipj)/2.0d0)
 527 1869c1505
 528 <             evdw=evdw+evdwij*sss
 529 ---
 530 >             evdw=evdw+evdwij
 531 1871,1872c1507,1508
 532 <             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 533 <             epsi=bb**2/aa
 534 ---
 535 >             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 536 >             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 537 1889,1891d1524
 538 < c            print '(2i4,6f8.4)',i,j,sss,sssgrad*
 539 < c     &      evdwij,fac,sigma(itypi,itypj),expon
 540 <             fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
 541 1894,1901d1526
 542 <             gg_lipi(3)=eps1*(eps2rt*eps2rt)
 543 <      &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
 544 <      & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
 545 <      &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
 546 <             gg_lipj(3)=ssgradlipj*gg_lipi(3)
 547 <             gg_lipi(3)=gg_lipi(3)*ssgradlipi
 548 < C            gg_lipi(3)=0.0d0
 549 < C            gg_lipj(3)=0.0d0
 550 1907,1908d1531
 551 <             endif
 552 <             ENDIF    ! dyn_ss
 553 1912,1914d1534
 554 < C      enddo          ! zshift
 555 < C      enddo          ! yshift
 556 < C      enddo          ! xshift
 557 1951,1985d1570
 558 <           xi=mod(xi,boxxsize)
 559 <           if (xi.lt.0) xi=xi+boxxsize
 560 <           yi=mod(yi,boxysize)
 561 <           if (yi.lt.0) yi=yi+boxysize
 562 <           zi=mod(zi,boxzsize)
 563 <           if (zi.lt.0) zi=zi+boxzsize
 564 < C define scaling factor for lipids
 565 <
 566 < C        if (positi.le.0) positi=positi+boxzsize
 567 < C        print *,i
 568 < C first for peptide groups
 569 < c for each residue check if it is in lipid or lipid water border area
 570 <        if ((zi.gt.bordlipbot)
 571 <      &.and.(zi.lt.bordliptop)) then
 572 < C the energy transfer exist
 573 <         if (zi.lt.buflipbot) then
 574 < C what fraction I am in
 575 <          fracinbuf=1.0d0-
 576 <      &        ((zi-bordlipbot)/lipbufthick)
 577 < C lipbufthick is thickenes of lipid buffore
 578 <          sslipi=sscalelip(fracinbuf)
 579 <          ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
 580 <         elseif (zi.gt.bufliptop) then
 581 <          fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
 582 <          sslipi=sscalelip(fracinbuf)
 583 <          ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
 584 <         else
 585 <          sslipi=1.0d0
 586 <          ssgradlipi=0.0
 587 <         endif
 588 <        else
 589 <          sslipi=0.0d0
 590 <          ssgradlipi=0.0
 591 <        endif
 592 <
 593 2022,2089c1607,1609
 594 < C            xj=c(1,nres+j)-xi
 595 < C            yj=c(2,nres+j)-yi
 596 < C            zj=c(3,nres+j)-zi
 597 <           xj=mod(xj,boxxsize)
 598 <           if (xj.lt.0) xj=xj+boxxsize
 599 <           yj=mod(yj,boxysize)
 600 <           if (yj.lt.0) yj=yj+boxysize
 601 <           zj=mod(zj,boxzsize)
 602 <           if (zj.lt.0) zj=zj+boxzsize
 603 <        if ((zj.gt.bordlipbot)
 604 <      &.and.(zj.lt.bordliptop)) then
 605 < C the energy transfer exist
 606 <         if (zj.lt.buflipbot) then
 607 < C what fraction I am in
 608 <          fracinbuf=1.0d0-
 609 <      &        ((zj-bordlipbot)/lipbufthick)
 610 < C lipbufthick is thickenes of lipid buffore
 611 <          sslipj=sscalelip(fracinbuf)
 612 <          ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
 613 <         elseif (zj.gt.bufliptop) then
 614 <          fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
 615 <          sslipj=sscalelip(fracinbuf)
 616 <          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
 617 <         else
 618 <          sslipj=1.0d0
 619 <          ssgradlipj=0.0
 620 <         endif
 621 <        else
 622 <          sslipj=0.0d0
 623 <          ssgradlipj=0.0
 624 <        endif
 625 <       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 626 <      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 627 <       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
 628 <      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 629 < C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5')
 630 < C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
 631 <       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 632 <       xj_safe=xj
 633 <       yj_safe=yj
 634 <       zj_safe=zj
 635 <       subchap=0
 636 <       do xshift=-1,1
 637 <       do yshift=-1,1
 638 <       do zshift=-1,1
 639 <           xj=xj_safe+xshift*boxxsize
 640 <           yj=yj_safe+yshift*boxysize
 641 <           zj=zj_safe+zshift*boxzsize
 642 <           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 643 <           if(dist_temp.lt.dist_init) then
 644 <             dist_init=dist_temp
 645 <             xj_temp=xj
 646 <             yj_temp=yj
 647 <             zj_temp=zj
 648 <             subchap=1
 649 <           endif
 650 <        enddo
 651 <        enddo
 652 <        enddo
 653 <        if (subchap.eq.1) then
 654 <           xj=xj_temp-xi
 655 <           yj=yj_temp-yi
 656 <           zj=zj_temp-zi
 657 <        else
 658 <           xj=xj_safe-xi
 659 <           yj=yj_safe-yi
 660 <           zj=zj_safe-zi
 661 <        endif
 662 ---
 663 >             xj=c(1,nres+j)-xi
 664 >             yj=c(2,nres+j)-yi
 665 >             zj=c(3,nres+j)-zi
 666 2110,2111c1630,1631
 667 <             e1=fac*fac*aa
 668 <             e2=fac*bb
 669 ---
 670 >             e1=fac*fac*aa(itypi,itypj)
 671 >             e2=fac*bb(itypi,itypj)
 672 2120,2121c1640,1641
 673 <             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
 674 <             epsi=bb**2/aa
 675 ---
 676 >             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 677 >             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 678 2135d1654
 679 <             fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
 680 2223d1741
 681 < cc      print *,'sss=',sss
 682 2242c1760
 683 <         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss
 684 ---
 685 >         gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
 686 2246c1764
 687 <         gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
 688 ---
 689 >         gvdwx(k,i)=gvdwx(k,i)-gg(k)
 690 2248,2249c1766,1767
 691 <      &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss
 692 <         gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
 693 ---
 694 >      &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 695 >         gvdwx(k,j)=gvdwx(k,j)+gg(k)
 696 2251c1769
 697 <      &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss
 698 ---
 699 >      &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
 700 2266,2267c1784,1785
 701 <         gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
 702 <         gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
 703 ---
 704 >         gvdwc(l,i)=gvdwc(l,i)-gg(l)
 705 >         gvdwc(l,j)=gvdwc(l,j)+gg(l)
 706 2369c1887
 707 < C      write(iout,*) 'In EELEC_soft_sphere'
 708 ---
 709 > cd      write(iout,*) 'In EELEC_soft_sphere'
 710 2384,2389d1901
 711 <           xmedi=mod(xmedi,boxxsize)
 712 <           if (xmedi.lt.0) xmedi=xmedi+boxxsize
 713 <           ymedi=mod(ymedi,boxysize)
 714 <           if (ymedi.lt.0) ymedi=ymedi+boxysize
 715 <           zmedi=mod(zmedi,boxzsize)
 716 <           if (zmedi.lt.0) zmedi=zmedi+boxzsize
 717 2403,2442c1915,1917
 718 <           xj=c(1,j)+0.5D0*dxj
 719 <           yj=c(2,j)+0.5D0*dyj
 720 <           zj=c(3,j)+0.5D0*dzj
 721 <           xj=mod(xj,boxxsize)
 722 <           if (xj.lt.0) xj=xj+boxxsize
 723 <           yj=mod(yj,boxysize)
 724 <           if (yj.lt.0) yj=yj+boxysize
 725 <           zj=mod(zj,boxzsize)
 726 <           if (zj.lt.0) zj=zj+boxzsize
 727 <       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
 728 <       xj_safe=xj
 729 <       yj_safe=yj
 730 <       zj_safe=zj
 731 <       isubchap=0
 732 <       do xshift=-1,1
 733 <       do yshift=-1,1
 734 <       do zshift=-1,1
 735 <           xj=xj_safe+xshift*boxxsize
 736 <           yj=yj_safe+yshift*boxysize
 737 <           zj=zj_safe+zshift*boxzsize
 738 <           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 739 <           if(dist_temp.lt.dist_init) then
 740 <             dist_init=dist_temp
 741 <             xj_temp=xj
 742 <             yj_temp=yj
 743 <             zj_temp=zj
 744 <             isubchap=1
 745 <           endif
 746 <        enddo
 747 <        enddo
 748 <        enddo
 749 <        if (isubchap.eq.1) then
 750 <           xj=xj_temp-xmedi
 751 <           yj=yj_temp-ymedi
 752 <           zj=zj_temp-zmedi
 753 <        else
 754 <           xj=xj_safe-xmedi
 755 <           yj=yj_safe-ymedi
 756 <           zj=zj_safe-zmedi
 757 <        endif
 758 ---
 759 >           xj=c(1,j)+0.5D0*dxj-xmedi
 760 >           yj=c(2,j)+0.5D0*dyj-ymedi
 761 >           zj=c(3,j)+0.5D0*dzj-zmedi
 762 2444,2445d1918
 763 <             sss=sscale(sqrt(rij))
 764 <             sssgrad=sscagrad(sqrt(rij))
 765 2453c1926
 766 <           evdw1=evdw1+evdw1ij*sss
 767 ---
 768 >           evdw1=evdw1+evdw1ij
 769 2457,2459c1930,1932
 770 <           ggg(1)=fac*xj*sssgrad
 771 <           ggg(2)=fac*yj*sssgrad
 772 <           ggg(3)=fac*zj*sssgrad
 773 ---
 774 >           ggg(1)=fac*xj
 775 >           ggg(2)=fac*yj
 776 >           ggg(3)=fac*zj
 777 2803c2276
 778 < c     &           +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
 779 ---
 780 > c     &           +bnew1(3,1,iti)*dsin(alpha(i))*cos(beta(i))
 781 2808c2281
 782 < c     &           +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
 783 ---
 784 > c     &           +bnew2(3,1,iti)*dsin(alpha(i))*dcos(beta(i))
 785 2815,2817c2288,2290
 786 < c     &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
 787 < c     &bnew1(2,1,iti)*cos(theta(i)),
 788 < c     &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0
 789 ---
 790 > c     &bnew1(1,1,iti)*dcos(theta(i)/2.0d0)/2.0d0,
 791 > c     &bnew1(2,1,iti)*dcos(theta(i)),
 792 > c     &bnew1(3,1,iti)*dsin(theta(i)/2.0d0)/2.0d0
 793 2834,2835c2307,2308
 794 < c        b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i))
 795 < c        b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
 796 ---
 797 > c        b1(2,iti)=bnew1(1,2,iti)*dsin(alpha(i))*dsin(beta(i))
 798 > c        b2(2,iti)=bnew2(1,2,iti)*dsin(alpha(i))*dsin(beta(i))
 799 2841d2313
 800 < c       write(iout,*)  'b1=',b1(1,i-2)
 801 2843a2316,2317
 802 > #ifdef PARMAT
 803 >       do i=ivec_start+2,ivec_end+2
 804 2844a2319,2321
 805 >       do i=3,nres+1
 806 > #endif
 807 > #endif
 808 2856,2874d2332
 809 <         b1(1,i-2)=b(3,iti)
 810 <         b1(2,i-2)=b(5,iti)
 811 <         b2(1,i-2)=b(2,iti)
 812 <         b2(2,i-2)=b(4,iti)
 813 <        b1tilde(1,i-2)=b1(1,i-2)
 814 <        b1tilde(2,i-2)=-b1(2,i-2)
 815 <        b2tilde(1,i-2)=b2(1,i-2)
 816 <        b2tilde(2,i-2)=-b2(2,i-2)
 817 <         EE(1,2,i-2)=eeold(1,2,iti)
 818 <         EE(2,1,i-2)=eeold(2,1,iti)
 819 <         EE(2,2,i-2)=eeold(2,2,iti)
 820 <         EE(1,1,i-2)=eeold(1,1,iti)
 821 <       enddo
 822 < #endif
 823 < #ifdef PARMAT
 824 <       do i=ivec_start+2,ivec_end+2
 825 < #else
 826 <       do i=3,nres+1
 827 < #endif
 828 2964,2965c2422,2423
 829 < c          write(iout,*) "co jest kurwa", iti, EE(1,1,i),EE(2,1,i),
 830 < c     &    EE(1,2,iti),EE(2,2,i)
 831 ---
 832 > c          write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti),
 833 > c     &    EE(1,2,iti),EE(2,2,iti)
 834 3010,3011c2468,2476
 835 < C        write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
 836 < c        write (iout,*) 'mu ',mu(:,i-2),i-2
 837 ---
 838 > #ifdef MUOUT
 839 >         write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1),
 840 >      &     rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2),
 841 >      &       -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2),
 842 >      &       dsqrt(b2(1,i-1)**2+b2(2,i-1)**2)
 843 >      &      +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2),
 844 >      &      ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itype2loc(iti))
 845 > #endif
 846 > cd        write (iout,*) 'mu ',mu(:,i-2)
 847 3298c2763
 848 < cd        iti = itortyp(itype(i))
 849 ---
 850 > cd        iti = itype2loc(itype(i))
 851 3335d2799
 852 <       include 'COMMON.SPLITELE'
 853 3418d2881
 854 < C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
 855 3420,3421d2882
 856 <         if (i.le.1) cycle
 857 < C        write(iout,*) "tu jest i",i
 858 3423,3434c2884
 859 < C changes suggested by Ana to avoid out of bounds
 860 <      & .or.((i+4).gt.nres)
 861 <      & .or.((i-1).le.0)
 862 < C end of changes by Ana
 863 <      &  .or. itype(i+2).eq.ntyp1
 864 <      &  .or. itype(i+3).eq.ntyp1) cycle
 865 <         if(i.gt.1)then
 866 <           if(itype(i-1).eq.ntyp1)cycle
 867 <         end if
 868 <         if(i.LT.nres-3)then
 869 <           if (itype(i+4).eq.ntyp1) cycle
 870 <         end if
 871 ---
 872 >      &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
 873 3444,3449d2893
 874 <           xmedi=mod(xmedi,boxxsize)
 875 <           if (xmedi.lt.0) xmedi=xmedi+boxxsize
 876 <           ymedi=mod(ymedi,boxysize)
 877 <           if (ymedi.lt.0) ymedi=ymedi+boxysize
 878 <           zmedi=mod(zmedi,boxzsize)
 879 <           if (zmedi.lt.0) zmedi=zmedi+boxzsize
 880 3456d2899
 881 <         if (i.le.1) cycle
 882 3458,3461d2900
 883 < C changes suggested by Ana to avoid out of bounds
 884 <      & .or.((i+5).gt.nres)
 885 <      & .or.((i-1).le.0)
 886 < C end of changes suggested by Ana
 887 3463,3467c2902
 888 <      &    .or. itype(i+4).eq.ntyp1
 889 <      &    .or. itype(i+5).eq.ntyp1
 890 <      &    .or. itype(i).eq.ntyp1
 891 <      &    .or. itype(i-1).eq.ntyp1
 892 <      &                             ) cycle
 893 ---
 894 >      &    .or. itype(i+4).eq.ntyp1) cycle
 895 3477,3508d2911
 896 < C Return atom into box, boxxsize is size of box in x dimension
 897 < c  194   continue
 898 < c        if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
 899 < c        if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
 900 < C Condition for being inside the proper box
 901 < c        if ((xmedi.gt.((0.5d0)*boxxsize)).or.
 902 < c     &       (xmedi.lt.((-0.5d0)*boxxsize))) then
 903 < c        go to 194
 904 < c        endif
 905 < c  195   continue
 906 < c        if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
 907 < c        if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
 908 < C Condition for being inside the proper box
 909 < c        if ((ymedi.gt.((0.5d0)*boxysize)).or.
 910 < c     &       (ymedi.lt.((-0.5d0)*boxysize))) then
 911 < c        go to 195
 912 < c        endif
 913 < c  196   continue
 914 < c        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
 915 < c        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
 916 < C Condition for being inside the proper box
 917 < c        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
 918 < c     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
 919 < c        go to 196
 920 < c        endif
 921 <           xmedi=mod(xmedi,boxxsize)
 922 <           if (xmedi.lt.0) xmedi=xmedi+boxxsize
 923 <           ymedi=mod(ymedi,boxysize)
 924 <           if (ymedi.lt.0) ymedi=ymedi+boxysize
 925 <           zmedi=mod(zmedi,boxzsize)
 926 <           if (zmedi.lt.0) zmedi=zmedi+boxzsize
 927 <
 928 3516,3519d2918
 929 < C Loop over all neighbouring boxes
 930 < C      do xshift=-1,1
 931 < C      do yshift=-1,1
 932 < C      do zshift=-1,1
 933 3523d2921
 934 < CTU KURWA
 935 3525,3534c2923,2924
 936 < C        do i=75,75
 937 <         if (i.le.1) cycle
 938 <         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 939 < C changes suggested by Ana to avoid out of bounds
 940 <      & .or.((i+2).gt.nres)
 941 <      & .or.((i-1).le.0)
 942 < C end of changes by Ana
 943 <      &  .or. itype(i+2).eq.ntyp1
 944 <      &  .or. itype(i-1).eq.ntyp1
 945 <      &                ) cycle
 946 ---
 947 > c       do i=7,7
 948 >         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
 949 3544,3579d2933
 950 <           xmedi=mod(xmedi,boxxsize)
 951 <           if (xmedi.lt.0) xmedi=xmedi+boxxsize
 952 <           ymedi=mod(ymedi,boxysize)
 953 <           if (ymedi.lt.0) ymedi=ymedi+boxysize
 954 <           zmedi=mod(zmedi,boxzsize)
 955 <           if (zmedi.lt.0) zmedi=zmedi+boxzsize
 956 < C          xmedi=xmedi+xshift*boxxsize
 957 < C          ymedi=ymedi+yshift*boxysize
 958 < C          zmedi=zmedi+zshift*boxzsize
 959 <
 960 < C Return tom into box, boxxsize is size of box in x dimension
 961 < c  164   continue
 962 < c        if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
 963 < c        if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
 964 < C Condition for being inside the proper box
 965 < c        if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
 966 < c     &       (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
 967 < c        go to 164
 968 < c        endif
 969 < c  165   continue
 970 < c        if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
 971 < c        if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
 972 < C Condition for being inside the proper box
 973 < c        if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
 974 < c     &       (ymedi.lt.((yshift-0.5d0)*boxysize))) then
 975 < c        go to 165
 976 < c        endif
 977 < c  166   continue
 978 < c        if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
 979 < c        if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
 980 < cC Condition for being inside the proper box
 981 < c        if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
 982 < c     &       (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
 983 < c        go to 166
 984 < c        endif
 985 <
 986 3582d2935
 987 < C I TU KURWA
 988 3584,3594c2937,2939
 989 < C          do j=16,17
 990 < C          write (iout,*) i,j
 991 <          if (j.le.1) cycle
 992 <           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
 993 < C changes suggested by Ana to avoid out of bounds
 994 <      & .or.((j+2).gt.nres)
 995 <      & .or.((j-1).le.0)
 996 < C end of changes by Ana
 997 <      & .or.itype(j+2).eq.ntyp1
 998 <      & .or.itype(j-1).eq.ntyp1
 999 <      &) cycle
1000 ---
1001 > c         do j=13,13
1002 > c          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
1003 >           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
1004 3599,3602d2943
1005 < C     enddo   ! zshift
1006 < C      enddo   ! yshift
1007 < C      enddo   ! xshift
1008 <
1009 3633,3634d2973
1010 <       include 'COMMON.SPLITELE'
1011 <       include 'COMMON.SHIELD'
1012 3670,3742c3009,3011
1013 < C          xj=c(1,j)+0.5D0*dxj-xmedi
1014 < C          yj=c(2,j)+0.5D0*dyj-ymedi
1015 < C          zj=c(3,j)+0.5D0*dzj-zmedi
1016 <           xj=c(1,j)+0.5D0*dxj
1017 <           yj=c(2,j)+0.5D0*dyj
1018 <           zj=c(3,j)+0.5D0*dzj
1019 <           xj=mod(xj,boxxsize)
1020 <           if (xj.lt.0) xj=xj+boxxsize
1021 <           yj=mod(yj,boxysize)
1022 <           if (yj.lt.0) yj=yj+boxysize
1023 <           zj=mod(zj,boxzsize)
1024 <           if (zj.lt.0) zj=zj+boxzsize
1025 <           if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
1026 <       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1027 <       xj_safe=xj
1028 <       yj_safe=yj
1029 <       zj_safe=zj
1030 <       isubchap=0
1031 <       do xshift=-1,1
1032 <       do yshift=-1,1
1033 <       do zshift=-1,1
1034 <           xj=xj_safe+xshift*boxxsize
1035 <           yj=yj_safe+yshift*boxysize
1036 <           zj=zj_safe+zshift*boxzsize
1037 <           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
1038 <           if(dist_temp.lt.dist_init) then
1039 <             dist_init=dist_temp
1040 <             xj_temp=xj
1041 <             yj_temp=yj
1042 <             zj_temp=zj
1043 <             isubchap=1
1044 <           endif
1045 <        enddo
1046 <        enddo
1047 <        enddo
1048 <        if (isubchap.eq.1) then
1049 <           xj=xj_temp-xmedi
1050 <           yj=yj_temp-ymedi
1051 <           zj=zj_temp-zmedi
1052 <        else
1053 <           xj=xj_safe-xmedi
1054 <           yj=yj_safe-ymedi
1055 <           zj=zj_safe-zmedi
1056 <        endif
1057 < C        if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
1058 < c  174   continue
1059 < c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
1060 < c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
1061 < C Condition for being inside the proper box
1062 < c        if ((xj.gt.((0.5d0)*boxxsize)).or.
1063 < c     &       (xj.lt.((-0.5d0)*boxxsize))) then
1064 < c        go to 174
1065 < c        endif
1066 < c  175   continue
1067 < c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
1068 < c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
1069 < C Condition for being inside the proper box
1070 < c        if ((yj.gt.((0.5d0)*boxysize)).or.
1071 < c     &       (yj.lt.((-0.5d0)*boxysize))) then
1072 < c        go to 175
1073 < c        endif
1074 < c  176   continue
1075 < c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
1076 < c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
1077 < C Condition for being inside the proper box
1078 < c        if ((zj.gt.((0.5d0)*boxzsize)).or.
1079 < c     &       (zj.lt.((-0.5d0)*boxzsize))) then
1080 < c        go to 176
1081 < c        endif
1082 < C        endif !endPBC condintion
1083 < C        xj=xj-xmedi
1084 < C        yj=yj-ymedi
1085 < C        zj=zj-zmedi
1086 ---
1087 >           xj=c(1,j)+0.5D0*dxj-xmedi
1088 >           yj=c(2,j)+0.5D0*dyj-ymedi
1089 >           zj=c(3,j)+0.5D0*dzj-zmedi
1090 3744,3747d3012
1091 <
1092 <             sss=sscale(sqrt(rij))
1093 <             sssgrad=sscagrad(sqrt(rij))
1094 < c            if (sss.gt.0.0d0) then
1095 3763c3028
1096 <           evdwij=(ev1+ev2)
1097 ---
1098 >           evdwij=ev1+ev2
1099 3766,3767c3031
1100 < C MARYSIA
1101 < C          eesij=(el1+el2)
1102 ---
1103 >           eesij=el1+el2
1104 3770,3780d3033
1105 <           if (shield_mode.gt.0) then
1106 < C          fac_shield(i)=0.4
1107 < C          fac_shield(j)=0.6
1108 <           el1=el1*fac_shield(i)**2*fac_shield(j)**2
1109 <           el2=el2*fac_shield(i)**2*fac_shield(j)**2
1110 <           eesij=(el1+el2)
1111 <           ees=ees+eesij
1112 <           else
1113 <           fac_shield(i)=1.0
1114 <           fac_shield(j)=1.0
1115 <           eesij=(el1+el2)
1116 3782,3783c3035
1117 <           endif
1118 <           evdw1=evdw1+evdwij*sss
1119 ---
1120 >           evdw1=evdw1+evdwij
1121 3793,3794c3045
1122 <               write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
1123 <      &fac_shield(i),fac_shield(j)
1124 ---
1125 >               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1126 3801c3052
1127 <           facvdw=-6*rrmij*(ev1+evdwij)*sss
1128 ---
1129 >           facvdw=-6*rrmij*(ev1+evdwij)
1130 3807d3057
1131 <
1132 3814,3882d3063
1133 <           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1134 <      &  (shield_mode.gt.0)) then
1135 < C          print *,i,j
1136 <           do ilist=1,ishield_list(i)
1137 <            iresshield=shield_list(ilist,i)
1138 <            do k=1,3
1139 <            rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1140 <      &      *2.0
1141 <            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1142 <      &              rlocshield
1143 <      & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
1144 <             gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1145 < C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1146 < C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1147 < C             if (iresshield.gt.i) then
1148 < C               do ishi=i+1,iresshield-1
1149 < C                gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1150 < C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1151 < C
1152 < C              enddo
1153 < C             else
1154 < C               do ishi=iresshield,i
1155 < C                gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1156 < C     & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
1157 < C
1158 < C               enddo
1159 < C              endif
1160 <            enddo
1161 <           enddo
1162 <           do ilist=1,ishield_list(j)
1163 <            iresshield=shield_list(ilist,j)
1164 <            do k=1,3
1165 <            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1166 <      &     *2.0
1167 <            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1168 <      &              rlocshield
1169 <      & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1170 <            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
1171 <
1172 < C     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
1173 < C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
1174 < C     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
1175 < C             if (iresshield.gt.j) then
1176 < C               do ishi=j+1,iresshield-1
1177 < C                gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
1178 < C     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
1179 < C
1180 < C               enddo
1181 < C            else
1182 < C               do ishi=iresshield,j
1183 < C                gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
1184 < C     & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
1185 < C               enddo
1186 < C              endif
1187 <            enddo
1188 <           enddo
1189 <
1190 <           do k=1,3
1191 <             gshieldc(k,i)=gshieldc(k,i)+
1192 <      &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
1193 <             gshieldc(k,j)=gshieldc(k,j)+
1194 <      &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
1195 <             gshieldc(k,i-1)=gshieldc(k,i-1)+
1196 <      &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
1197 <             gshieldc(k,j-1)=gshieldc(k,j-1)+
1198 <      &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
1199 <
1200 <            enddo
1201 <            endif
1202 3889d3069
1203 < C           print *,"before", gelc_long(1,i), gelc_long(1,j)
1204 3892d3071
1205 < C     &                    +grad_shield(k,j)*eesij/fac_shield(j)
1206 3894,3898d3072
1207 < C     &                    +grad_shield(k,i)*eesij/fac_shield(i)
1208 < C            gelc_long(k,i-1)=gelc_long(k,i-1)
1209 < C     &                    +grad_shield(k,i)*eesij/fac_shield(i)
1210 < C            gelc_long(k,j-1)=gelc_long(k,j-1)
1211 < C     &                    +grad_shield(k,j)*eesij/fac_shield(j)
1212 3900,3901d3073
1213 < C           print *,"bafter", gelc_long(1,i), gelc_long(1,j)
1214 <
1215 3910,3918c3082,3084
1216 <           if (sss.gt.0.0) then
1217 <           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
1218 <           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
1219 <           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
1220 <           else
1221 <           ggg(1)=0.0
1222 <           ggg(2)=0.0
1223 <           ggg(3)=0.0
1224 <           endif
1225 ---
1226 >           ggg(1)=facvdw*xj
1227 >           ggg(2)=facvdw*yj
1228 >           ggg(3)=facvdw*zj
1229 3938,3940c3104,3105
1230 < C MARYSIA
1231 <           facvdw=(ev1+evdwij)*sss
1232 <           facel=(el1+eesij)
1233 ---
1234 >           facvdw=ev1+evdwij
1235 >           facel=el1+eesij
1236 3950d3114
1237 < C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j)
1238 3952d3115
1239 < C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j)
1240 3954d3116
1241 < C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j)
1242 3974,3976c3136,3138
1243 <           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
1244 <           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
1245 <           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
1246 ---
1247 >           ggg(1)=facvdw*xj
1248 >           ggg(2)=facvdw*yj
1249 >           ggg(3)=facvdw*zj
1250 3997,3998c3159
1251 <             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1252 <      &      fac_shield(i)**2*fac_shield(j)**2
1253 ---
1254 >             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1255 4014d3174
1256 < C                     print *,"before22", gelc_long(1,i), gelc_long(1,j)
1257 4017,4019c3177,3178
1258 <      &           +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1259 <      &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
1260 <      &           *fac_shield(i)**2*fac_shield(j)**2
1261 ---
1262 >      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1263 >      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
1264 4021,4023c3180,3181
1265 <      &           +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1266 <      &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
1267 <      &           *fac_shield(i)**2*fac_shield(j)**2
1268 ---
1269 >      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1270 >      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1271 4027,4030d3184
1272 < C           print *,"before33", gelc_long(1,i), gelc_long(1,j)
1273 <
1274 < C MARYSIA
1275 < c          endif !sscale
1276 4055d3208
1277 < c              write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
1278 4058c3211
1279 < c             write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
1280 ---
1281 > c             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
1282 4061c3214
1283 < c             write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
1284 ---
1285 > c             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
1286 4083c3236
1287 < cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1288 ---
1289 > cd     &     i,itype2loc(itype(i)),j,itype2loc(itype(j)),a22,a23,a32,a33
1290 4225a3379
1291 > c           if ((i.eq.8).and.(j.eq.14)) then
1292 4229,4286d3382
1293 <           if (shield_mode.eq.0) then
1294 <            fac_shield(i)=1.0
1295 <            fac_shield(j)=1.0
1296 < C          else
1297 < C           fac_shield(i)=0.4
1298 < C           fac_shield(j)=0.6
1299 <           endif
1300 <           eel_loc_ij=eel_loc_ij
1301 <      &    *fac_shield(i)*fac_shield(j)
1302 < C Now derivative over eel_loc
1303 <           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1304 <      &  (shield_mode.gt.0)) then
1305 < C          print *,i,j
1306 <
1307 <           do ilist=1,ishield_list(i)
1308 <            iresshield=shield_list(ilist,i)
1309 <            do k=1,3
1310 <            rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
1311 <      &                                          /fac_shield(i)
1312 < C     &      *2.0
1313 <            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
1314 <      &              rlocshield
1315 <      & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
1316 <             gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
1317 <      &      +rlocshield
1318 <            enddo
1319 <           enddo
1320 <           do ilist=1,ishield_list(j)
1321 <            iresshield=shield_list(ilist,j)
1322 <            do k=1,3
1323 <            rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
1324 <      &                                       /fac_shield(j)
1325 < C     &     *2.0
1326 <            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
1327 <      &              rlocshield
1328 <      & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
1329 <            gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
1330 <      &             +rlocshield
1331 <
1332 <            enddo
1333 <           enddo
1334 <
1335 <           do k=1,3
1336 <             gshieldc_ll(k,i)=gshieldc_ll(k,i)+
1337 <      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
1338 <             gshieldc_ll(k,j)=gshieldc_ll(k,j)+
1339 <      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
1340 <             gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
1341 <      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
1342 <             gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
1343 <      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
1344 <            enddo
1345 <            endif
1346 <
1347 <
1348 < c          write (iout,*) 'i',i,' j',j,itype(i),itype(j),
1349 < c     &                     ' eel_loc_ij',eel_loc_ij
1350 < C          write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
1351 4289c3385
1352 <          geel_loc_ij=(a22*gmuij1(1)
1353 ---
1354 >          geel_loc_ij=a22*gmuij1(1)
1355 4292,4293c3388
1356 <      &     +a33*gmuij1(4))
1357 <      &    *fac_shield(i)*fac_shield(j)
1358 ---
1359 >      &     +a33*gmuij1(4)
1360 4302,4305c3397
1361 <          geel_loc_ij=
1362 <      &     a22*gmuij2(1)
1363 <      &     +a23*gmuij2(2)
1364 <      &     +a32*gmuij2(3)
1365 ---
1366 >          geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
1367 4309,4314c3401
1368 <      &    *fac_shield(i)*fac_shield(j)
1369 <
1370 < c  Derivative over j residue
1371 <          geel_loc_ji=a22*gmuji1(1)
1372 <      &     +a23*gmuji1(2)
1373 <      &     +a32*gmuji1(3)
1374 ---
1375 >          geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
1376 4320c3407
1377 <         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
1378 ---
1379 >          gloc(nphi+j,icg)=gloc(nphi+j,icg)+
1380 4322,4327c3409
1381 <      &    *fac_shield(i)*fac_shield(j)
1382 <
1383 <          geel_loc_ji=
1384 <      &     +a22*gmuji2(1)
1385 <      &     +a23*gmuji2(2)
1386 <      &     +a32*gmuji2(3)
1387 ---
1388 >          geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
1389 4334d3415
1390 <      &    *fac_shield(i)*fac_shield(j)
1391 4340,4342c3421
1392 < c           if (eel_loc_ij.ne.0)
1393 < c     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
1394 < c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
1395 ---
1396 > c              write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
1397 4348,4351c3427,3428
1398 <      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1399 <      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
1400 <      &    *fac_shield(i)*fac_shield(j)
1401 <
1402 ---
1403 >      &            a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
1404 >      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
1405 4353,4355c3430,3431
1406 <      &           (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1407 <      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
1408 <      &    *fac_shield(i)*fac_shield(j)
1409 ---
1410 >      &            a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
1411 >      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
1412 4358,4360c3434,3435
1413 <             ggg(l)=(agg(l,1)*muij(1)+
1414 <      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
1415 <      &    *fac_shield(i)*fac_shield(j)
1416 ---
1417 >             ggg(l)=agg(l,1)*muij(1)+
1418 >      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
1419 4374,4389c3449,3456
1420 <             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
1421 <      &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
1422 <      &    *fac_shield(i)*fac_shield(j)
1423 <
1424 <             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
1425 <      &     aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
1426 <      &    *fac_shield(i)*fac_shield(j)
1427 <
1428 <             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
1429 <      &       aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
1430 <      &    *fac_shield(i)*fac_shield(j)
1431 <
1432 <             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
1433 <      &     aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
1434 <      &    *fac_shield(i)*fac_shield(j)
1435 <
1436 ---
1437 >             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
1438 >      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
1439 >             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
1440 >      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
1441 >             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
1442 >      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
1443 >             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
1444 >      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
1445 4390a3458
1446 > c          endif
1447 4468,4475d3535
1448 <                 if (shield_mode.eq.0) then
1449 <                 fac_shield(i)=1.0d0
1450 <                 fac_shield(j)=1.0d0
1451 <                 else
1452 <                 ees0plist(num_conti,i)=j
1453 < C                fac_shield(i)=0.4d0
1454 < C                fac_shield(j)=0.6d0
1455 <                 endif
1456 4477d3536
1457 <      &          *fac_shield(i)*fac_shield(j)
1458 4479d3537
1459 <      &          *fac_shield(i)*fac_shield(j)
1460 4547,4548d3604
1461 <      &          *fac_shield(i)*fac_shield(j)
1462 <
1463 4552,4553d3607
1464 <      &          *fac_shield(i)*fac_shield(j)
1465 <
1466 4555,4556d3608
1467 <      &          *fac_shield(i)*fac_shield(j)
1468 <
1469 4560,4561d3611
1470 <      &          *fac_shield(i)*fac_shield(j)
1471 <
1472 4565,4566d3614
1473 <      &          *fac_shield(i)*fac_shield(j)
1474 <
1475 4568,4569d3615
1476 <      &          *fac_shield(i)*fac_shield(j)
1477 <
1478 4621d3666
1479 <       include 'COMMON.SHIELD'
1480 4662,4668d3706
1481 <         if (shield_mode.eq.0) then
1482 <         fac_shield(i)=1.0
1483 <         fac_shield(j)=1.0
1484 < C        else
1485 < C        fac_shield(i)=0.4
1486 < C        fac_shield(j)=0.6
1487 <         endif
1488 4670,4672d3707
1489 <      &  *fac_shield(i)*fac_shield(j)
1490 <         eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
1491 <      &  *fac_shield(i)*fac_shield(j)
1492 4676d3710
1493 <      &   *fac_shield(i)*fac_shield(j)
1494 4679,4709d3712
1495 <      &   *fac_shield(i)*fac_shield(j)
1496 <
1497 <
1498 < C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1499 < C Derivatives in shield mode
1500 <           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1501 <      &  (shield_mode.gt.0)) then
1502 < C          print *,i,j
1503 <
1504 <           do ilist=1,ishield_list(i)
1505 <            iresshield=shield_list(ilist,i)
1506 <            do k=1,3
1507 <            rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i)
1508 < C     &      *2.0
1509 <            gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
1510 <      &              rlocshield
1511 <      & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i)
1512 <             gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
1513 <      &      +rlocshield
1514 <            enddo
1515 <           enddo
1516 <           do ilist=1,ishield_list(j)
1517 <            iresshield=shield_list(ilist,j)
1518 <            do k=1,3
1519 <            rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j)
1520 < C     &     *2.0
1521 <            gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
1522 <      &              rlocshield
1523 <      & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j)
1524 <            gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
1525 <      &             +rlocshield
1526 4711,4726c3714,3715
1527 <            enddo
1528 <           enddo
1529 <
1530 <           do k=1,3
1531 <             gshieldc_t3(k,i)=gshieldc_t3(k,i)+
1532 <      &              grad_shield(k,i)*eello_t3/fac_shield(i)
1533 <             gshieldc_t3(k,j)=gshieldc_t3(k,j)+
1534 <      &              grad_shield(k,j)*eello_t3/fac_shield(j)
1535 <             gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+
1536 <      &              grad_shield(k,i)*eello_t3/fac_shield(i)
1537 <             gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+
1538 <      &              grad_shield(k,j)*eello_t3/fac_shield(j)
1539 <            enddo
1540 <            endif
1541 <
1542 < C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1543 ---
1544 >         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1545 >      &          'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
1546 4735d3723
1547 <      &   *fac_shield(i)*fac_shield(j)
1548 4742d3729
1549 <      &   *fac_shield(i)*fac_shield(j)
1550 4756,4757d3742
1551 <      &   *fac_shield(i)*fac_shield(j)
1552 <
1553 4765d3749
1554 <      &   *fac_shield(i)*fac_shield(j)
1555 4773d3756
1556 <      &   *fac_shield(i)*fac_shield(j)
1557 4781d3763
1558 <      &   *fac_shield(i)*fac_shield(j)
1559 4802d3783
1560 <       include 'COMMON.SHIELD'
1561 4903,4961d3883
1562 <         if (shield_mode.eq.0) then
1563 <         fac_shield(i)=1.0
1564 <         fac_shield(j)=1.0
1565 < C        else
1566 < C        fac_shield(i)=0.6
1567 < C        fac_shield(j)=0.4
1568 <         endif
1569 <         eello_turn4=eello_turn4-(s1+s2+s3)
1570 <      &  *fac_shield(i)*fac_shield(j)
1571 <         eello_t4=-(s1+s2+s3)
1572 <      &  *fac_shield(i)*fac_shield(j)
1573 < c             write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
1574 <         if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)')
1575 <      &      'eturn4',i,j,-(s1+s2+s3),s1,s2,s3
1576 < C Now derivative over shield:
1577 <           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1578 <      &  (shield_mode.gt.0)) then
1579 < C          print *,i,j
1580 <
1581 <           do ilist=1,ishield_list(i)
1582 <            iresshield=shield_list(ilist,i)
1583 <            do k=1,3
1584 <            rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
1585 < C     &      *2.0
1586 <            gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
1587 <      &              rlocshield
1588 <      & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)
1589 <             gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
1590 <      &      +rlocshield
1591 <            enddo
1592 <           enddo
1593 <           do ilist=1,ishield_list(j)
1594 <            iresshield=shield_list(ilist,j)
1595 <            do k=1,3
1596 <            rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j)
1597 < C     &     *2.0
1598 <            gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
1599 <      &              rlocshield
1600 <      & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j)
1601 <            gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
1602 <      &             +rlocshield
1603 <
1604 <            enddo
1605 <           enddo
1606 <
1607 <           do k=1,3
1608 <             gshieldc_t4(k,i)=gshieldc_t4(k,i)+
1609 <      &              grad_shield(k,i)*eello_t4/fac_shield(i)
1610 <             gshieldc_t4(k,j)=gshieldc_t4(k,j)+
1611 <      &              grad_shield(k,j)*eello_t4/fac_shield(j)
1612 <             gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+
1613 <      &              grad_shield(k,i)*eello_t4/fac_shield(i)
1614 <             gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+
1615 <      &              grad_shield(k,j)*eello_t4/fac_shield(j)
1616 <            enddo
1617 <            endif
1618 <
1619 <
1620 <
1621 4963,4966c3885
1622 <
1623 <
1624 < cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
1625 < cd     &    ' eello_turn4_num',8*eello_turn4_num
1626 ---
1627 >         eello_turn4=eello_turn4-(s1+s2+s3)
1628 4970d3888
1629 <      &  *fac_shield(i)*fac_shield(j)
1630 4973,4974d3890
1631 <      &  *fac_shield(i)*fac_shield(j)
1632 <
1633 4977,4978d3892
1634 <      &  *fac_shield(i)*fac_shield(j)
1635 <
1636 4994d3907
1637 <      &  *fac_shield(i)*fac_shield(j)
1638 5003d3915
1639 <      &  *fac_shield(i)*fac_shield(j)
1640 5015d3926
1641 <      &  *fac_shield(i)*fac_shield(j)
1642 5035d3945
1643 <      &  *fac_shield(i)*fac_shield(j)
1644 5054d3963
1645 <      &  *fac_shield(i)*fac_shield(j)
1646 5069d3977
1647 <      &  *fac_shield(i)*fac_shield(j)
1648 5084d3991
1649 <      &  *fac_shield(i)*fac_shield(j)
1650 5100d4006
1651 <      &  *fac_shield(i)*fac_shield(j)
1652 5161,5163d4066
1653 < C      do xshift=-1,1
1654 < C      do yshift=-1,1
1655 < C      do zshift=-1,1
1656 5170,5203c4073
1657 < C Return atom into box, boxxsize is size of box in x dimension
1658 < c  134   continue
1659 < c        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
1660 < c        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
1661 < C Condition for being inside the proper box
1662 < c        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
1663 < c     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
1664 < c        go to 134
1665 < c        endif
1666 < c  135   continue
1667 < c        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
1668 < c        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
1669 < C Condition for being inside the proper box
1670 < c        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
1671 < c     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
1672 < c        go to 135
1673 < c c       endif
1674 < c  136   continue
1675 < c        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
1676 < c        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
1677 < cC Condition for being inside the proper box
1678 < c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
1679 < c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
1680 < c        go to 136
1681 < c        endif
1682 <           xi=mod(xi,boxxsize)
1683 <           if (xi.lt.0) xi=xi+boxxsize
1684 <           yi=mod(yi,boxysize)
1685 <           if (yi.lt.0) yi=yi+boxysize
1686 <           zi=mod(zi,boxzsize)
1687 <           if (zi.lt.0) zi=zi+boxzsize
1688 < C          xi=xi+xshift*boxxsize
1689 < C          yi=yi+yshift*boxysize
1690 < C          zi=zi+zshift*boxzsize
1691 ---
1692 >
1693 5214,5280c4084,4086
1694 <           xj=c(1,j)
1695 <           yj=c(2,j)
1696 <           zj=c(3,j)
1697 < c  174   continue
1698 < c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
1699 < c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
1700 < C Condition for being inside the proper box
1701 < c        if ((xj.gt.((0.5d0)*boxxsize)).or.
1702 < c     &       (xj.lt.((-0.5d0)*boxxsize))) then
1703 < c        go to 174
1704 < c        endif
1705 < c  175   continue
1706 < c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
1707 < c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
1708 < cC Condition for being inside the proper box
1709 < c        if ((yj.gt.((0.5d0)*boxysize)).or.
1710 < c     &       (yj.lt.((-0.5d0)*boxysize))) then
1711 < c        go to 175
1712 < c        endif
1713 < c  176   continue
1714 < c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
1715 < c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
1716 < C Condition for being inside the proper box
1717 < c        if ((zj.gt.((0.5d0)*boxzsize)).or.
1718 < c     &       (zj.lt.((-0.5d0)*boxzsize))) then
1719 < c        go to 176
1720 <           xj=mod(xj,boxxsize)
1721 <           if (xj.lt.0) xj=xj+boxxsize
1722 <           yj=mod(yj,boxysize)
1723 <           if (yj.lt.0) yj=yj+boxysize
1724 <           zj=mod(zj,boxzsize)
1725 <           if (zj.lt.0) zj=zj+boxzsize
1726 <       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1727 <       xj_safe=xj
1728 <       yj_safe=yj
1729 <       zj_safe=zj
1730 <       subchap=0
1731 <       do xshift=-1,1
1732 <       do yshift=-1,1
1733 <       do zshift=-1,1
1734 <           xj=xj_safe+xshift*boxxsize
1735 <           yj=yj_safe+yshift*boxysize
1736 <           zj=zj_safe+zshift*boxzsize
1737 <           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1738 <           if(dist_temp.lt.dist_init) then
1739 <             dist_init=dist_temp
1740 <             xj_temp=xj
1741 <             yj_temp=yj
1742 <             zj_temp=zj
1743 <             subchap=1
1744 <           endif
1745 <        enddo
1746 <        enddo
1747 <        enddo
1748 <        if (subchap.eq.1) then
1749 <           xj=xj_temp-xi
1750 <           yj=yj_temp-yi
1751 <           zj=zj_temp-zi
1752 <        else
1753 <           xj=xj_safe-xi
1754 <           yj=yj_safe-yi
1755 <           zj=zj_safe-zi
1756 <        endif
1757 < c c       endif
1758 < C          xj=xj-xi
1759 < C          yj=yj-yi
1760 < C          zj=zj-zi
1761 ---
1762 >           xj=c(1,j)-xi
1763 >           yj=c(2,j)-yi
1764 >           zj=c(3,j)-zi
1765 5282d4087
1766 <
1767 5333,5335d4137
1768 < C      enddo !zshift
1769 < C      enddo !yshift
1770 < C      enddo !xshift
1771 5356d4157
1772 <       include 'COMMON.SPLITELE'
1773 5360d4160
1774 < c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
1775 5363,5365d4162
1776 < C      do xshift=-1,1
1777 < C      do yshift=-1,1
1778 < C      do zshift=-1,1
1779 5372,5392d4168
1780 <           xi=mod(xi,boxxsize)
1781 <           if (xi.lt.0) xi=xi+boxxsize
1782 <           yi=mod(yi,boxysize)
1783 <           if (yi.lt.0) yi=yi+boxysize
1784 <           zi=mod(zi,boxzsize)
1785 <           if (zi.lt.0) zi=zi+boxzsize
1786 < c          xi=xi+xshift*boxxsize
1787 < c          yi=yi+yshift*boxysize
1788 < c          zi=zi+zshift*boxzsize
1789 < c        print *,xi,yi,zi,'polozenie i'
1790 < C Return atom into box, boxxsize is size of box in x dimension
1791 < c  134   continue
1792 < c        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
1793 < c        if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
1794 < C Condition for being inside the proper box
1795 < c        if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
1796 < c     &       (xi.lt.((xshift-0.5d0)*boxxsize))) then
1797 < c        go to 134
1798 < c        endif
1799 < c  135   continue
1800 < c          print *,xi,boxxsize,"pierwszy"
1801 5394,5408d4169
1802 < c        if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
1803 < c        if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
1804 < C Condition for being inside the proper box
1805 < c        if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
1806 < c     &       (yi.lt.((yshift-0.5d0)*boxysize))) then
1807 < c        go to 135
1808 < c        endif
1809 < c  136   continue
1810 < c        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
1811 < c        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
1812 < C Condition for being inside the proper box
1813 < c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
1814 < c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
1815 < c        go to 136
1816 < c        endif
1817 5419,5484c4180,4182
1818 <           xj=c(1,j)
1819 <           yj=c(2,j)
1820 <           zj=c(3,j)
1821 <           xj=mod(xj,boxxsize)
1822 <           if (xj.lt.0) xj=xj+boxxsize
1823 <           yj=mod(yj,boxysize)
1824 <           if (yj.lt.0) yj=yj+boxysize
1825 <           zj=mod(zj,boxzsize)
1826 <           if (zj.lt.0) zj=zj+boxzsize
1827 < c  174   continue
1828 < c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
1829 < c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
1830 < C Condition for being inside the proper box
1831 < c        if ((xj.gt.((0.5d0)*boxxsize)).or.
1832 < c     &       (xj.lt.((-0.5d0)*boxxsize))) then
1833 < c        go to 174
1834 < c        endif
1835 < c  175   continue
1836 < c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
1837 < c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
1838 < cC Condition for being inside the proper box
1839 < c        if ((yj.gt.((0.5d0)*boxysize)).or.
1840 < c     &       (yj.lt.((-0.5d0)*boxysize))) then
1841 < c        go to 175
1842 < c        endif
1843 < c  176   continue
1844 < c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
1845 < c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
1846 < C Condition for being inside the proper box
1847 < c        if ((zj.gt.((0.5d0)*boxzsize)).or.
1848 < c     &       (zj.lt.((-0.5d0)*boxzsize))) then
1849 < c        go to 176
1850 < c        endif
1851 < CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
1852 <       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1853 <       xj_safe=xj
1854 <       yj_safe=yj
1855 <       zj_safe=zj
1856 <       subchap=0
1857 <       do xshift=-1,1
1858 <       do yshift=-1,1
1859 <       do zshift=-1,1
1860 <           xj=xj_safe+xshift*boxxsize
1861 <           yj=yj_safe+yshift*boxysize
1862 <           zj=zj_safe+zshift*boxzsize
1863 <           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
1864 <           if(dist_temp.lt.dist_init) then
1865 <             dist_init=dist_temp
1866 <             xj_temp=xj
1867 <             yj_temp=yj
1868 <             zj_temp=zj
1869 <             subchap=1
1870 <           endif
1871 <        enddo
1872 <        enddo
1873 <        enddo
1874 <        if (subchap.eq.1) then
1875 <           xj=xj_temp-xi
1876 <           yj=yj_temp-yi
1877 <           zj=zj_temp-zi
1878 <        else
1879 <           xj=xj_safe-xi
1880 <           yj=yj_safe-yi
1881 <           zj=zj_safe-zi
1882 <        endif
1883 < c          print *,xj,yj,zj,'polozenie j'
1884 ---
1885 >           xj=c(1,j)-xi
1886 >           yj=c(2,j)-yi
1887 >           zj=c(3,j)-zi
1888 5486,5491d4183
1889 < c          print *,rrij
1890 <           sss=sscale(1.0d0/(dsqrt(rrij)))
1891 < c          print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz'
1892 < c          if (sss.eq.0) print *,'czasem jest OK'
1893 <           if (sss.le.0.0d0) cycle
1894 <           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
1895 5498c4190
1896 <             evdw2_14=evdw2_14+(e1+e2)*sss
1897 ---
1898 >             evdw2_14=evdw2_14+e1+e2
1899 5501c4193
1900 <           evdw2=evdw2+evdwij*sss
1901 ---
1902 >           evdw2=evdw2+evdwij
1903 5508,5509c4200
1904 <           fac=-(evdwij+e1)*rrij*sss
1905 <           fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
1906 ---
1907 >           fac=-(evdwij+e1)*rrij
1908 5544,5545c4235
1909 < c        endif !endif for sscale cutoff
1910 <         enddo ! j
1911 ---
1912 >         enddo
1913 5549,5551d4238
1914 < c      enddo !zshift
1915 < c      enddo !yshift
1916 < c      enddo !xshift
1917 5583d4269
1918 <       include 'COMMON.CONTROL'
1919 5586,5589d4271
1920 <       do i=1,3
1921 <        ggg(i)=0.0d0
1922 <       enddo
1923 < C      write (iout,*) ,"link_end",link_end,constr_dist
1924 5606,5607c4288
1925 < c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
1926 < c     &    dhpb(i),dhpb1(i),forcon(i)
1927 ---
1928 > cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
1929 5610,5616c4291
1930 < C        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
1931 < C     & iabs(itype(jjj)).eq.1) then
1932 < cmc        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
1933 < C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
1934 <         if (.not.dyn_ss .and. i.le.nss) then
1935 < C 15/02/13 CC dynamic SSbond - additional check
1936 <          if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
1937 ---
1938 >         if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
1939 5620d4294
1940 <          endif
1941 5622,5663d4295
1942 < cd   &   ' waga=',waga,' fac=',fac
1943 <         else if (ii.gt.nres .and. jj.gt.nres) then
1944 < c Restraints from contact prediction
1945 <           dd=dist(ii,jj)
1946 <           if (constr_dist.eq.11) then
1947 <             ehpb=ehpb+fordepth(i)**4.0d0
1948 <      &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
1949 <             fac=fordepth(i)**4.0d0
1950 <      &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
1951 <           if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
1952 <      &    ehpb,fordepth(i),dd
1953 <            else
1954 <           if (dhpb1(i).gt.0.0d0) then
1955 <             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
1956 <             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
1957 < c            write (iout,*) "beta nmr",
1958 < c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
1959 <           else
1960 <             dd=dist(ii,jj)
1961 <             rdis=dd-dhpb(i)
1962 < C Get the force constant corresponding to this distance.
1963 <             waga=forcon(i)
1964 < C Calculate the contribution to energy.
1965 <             ehpb=ehpb+waga*rdis*rdis
1966 < c            write (iout,*) "beta reg",dd,waga*rdis*rdis
1967 < C
1968 < C Evaluate gradient.
1969 < C
1970 <             fac=waga*rdis/dd
1971 <           endif
1972 <           endif
1973 <           do j=1,3
1974 <             ggg(j)=fac*(c(j,jj)-c(j,ii))
1975 <           enddo
1976 <           do j=1,3
1977 <             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
1978 <             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
1979 <           enddo
1980 <           do k=1,3
1981 <             ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
1982 <             ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
1983 <           enddo
1984 5667,5682c4299,4300
1985 <           dd=dist(ii,jj)
1986 <           if (constr_dist.eq.11) then
1987 <             ehpb=ehpb+fordepth(i)**4.0d0
1988 <      &           *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
1989 <             fac=fordepth(i)**4.0d0
1990 <      &           *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
1991 <           if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
1992 <      &    ehpb,fordepth(i),dd
1993 <            else
1994 <           if (dhpb1(i).gt.0.0d0) then
1995 <             ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
1996 <             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
1997 < c            write (iout,*) "alph nmr",
1998 < c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
1999 <           else
2000 <             rdis=dd-dhpb(i)
2001 ---
2002 >         dd=dist(ii,jj)
2003 >         rdis=dd-dhpb(i)
2004 5684c4302
2005 <             waga=forcon(i)
2006 ---
2007 >         waga=forcon(i)
2008 5686,5687c4304
2009 <             ehpb=ehpb+waga*rdis*rdis
2010 < c            write (iout,*) "alpha reg",dd,waga*rdis*rdis
2011 ---
2012 >         ehpb=ehpb+waga*rdis*rdis
2013 5691,5696c4308,4313
2014 <             fac=waga*rdis/dd
2015 <           endif
2016 <           endif
2017 <             do j=1,3
2018 <               ggg(j)=fac*(c(j,jj)-c(j,ii))
2019 <             enddo
2020 ---
2021 >         fac=waga*rdis/dd
2022 > cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
2023 > cd   &   ' waga=',waga,' fac=',fac
2024 >         do j=1,3
2025 >           ggg(j)=fac*(c(j,jj)-c(j,ii))
2026 >         enddo
2027 5700c4317
2028 <           if (iii.lt.ii) then
2029 ---
2030 >         if (iii.lt.ii) then
2031 5705c4322
2032 <           endif
2033 ---
2034 >         endif
2035 5711,5714c4328,4331
2036 <           do k=1,3
2037 <             ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
2038 <             ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
2039 <           enddo
2040 ---
2041 >         do k=1,3
2042 >           ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
2043 >           ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
2044 >         enddo
2045 5717c4334
2046 <       if (constr_dist.ne.11) ehpb=0.5D0*ehpb
2047 ---
2048 >       ehpb=0.5D0*ehpb
2049 5830,5844c4447,4455
2050 <         if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
2051 < c          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
2052 < c          do j=1,3
2053 < c          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
2054 < c     &      *dc(j,i-1)/vbld(i)
2055 < c          enddo
2056 < c          if (energy_dec) write(iout,*)
2057 < c     &       "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
2058 < c        else
2059 < C       Checking if it involves dummy (NH3+ or COO-) group
2060 <          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
2061 < C YES   vbldpDUM is the equlibrium length of spring for Dummy atom
2062 <         diff = vbld(i)-vbldpDUM
2063 <          else
2064 < C NO    vbldp0 is the equlibrium lenght of spring for peptide group
2065 ---
2066 >         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
2067 >           estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
2068 >           do j=1,3
2069 >           gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
2070 >      &      *dc(j,i-1)/vbld(i)
2071 >           enddo
2072 >           if (energy_dec) write(iout,*)
2073 >      &       "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
2074 >         else
2075 5846,5847c4457
2076 <          endif
2077 <         if (energy_dec)    write (iout,'(a7,i5,4f7.3)')
2078 ---
2079 >         if (energy_dec) write (iout,*)
2080 5854c4464
2081 < c        endif
2082 ---
2083 >         endif
2084 5866c4476
2085 <             if (energy_dec)  write (iout,*)
2086 ---
2087 >             if (energy_dec) write (iout,*)
2088 5908c4518
2089 <       subroutine ebend(etheta,ethetacnstr)
2090 ---
2091 >       subroutine ebend(etheta)
2092 5925d4534
2093 <       include 'COMMON.TORCNSTR'
2094 5936,5937c4545
2095 <         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
2096 <      &  .or.itype(i).eq.ntyp1) cycle
2097 ---
2098 >         if (itype(i-1).eq.ntyp1) cycle
2099 5954c4562
2100 <         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
2101 ---
2102 >         if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
2103 5967c4575
2104 <         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
2105 ---
2106 >         if (i.lt.nres .and. itype(i).ne.ntyp1) then
2107 5975d4582
2108 < #endif
2109 5976a4584
2110 > #endif
2111 5994d4601
2112 < c         write(iout,*) 'chuj tu', y(k),z(k)
2113 6031,6032c4638,4639
2114 <         if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
2115 <      &      'ebend',i,ethetai,theta(i),itype(i)
2116 ---
2117 >         if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
2118 >      &      'ebend',i,ethetai
2119 6037,6064d4643
2120 <       ethetacnstr=0.0d0
2121 < C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
2122 <       do i=ithetaconstr_start,ithetaconstr_end
2123 <         itheta=itheta_constr(i)
2124 <         thetiii=theta(itheta)
2125 <         difi=pinorm(thetiii-theta_constr0(i))
2126 <         if (difi.gt.theta_drange(i)) then
2127 <           difi=difi-theta_drange(i)
2128 <           ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
2129 <           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
2130 <      &    +for_thet_constr(i)*difi**3
2131 <         else if (difi.lt.-drange(i)) then
2132 <           difi=difi+drange(i)
2133 <           ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
2134 <           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
2135 <      &    +for_thet_constr(i)*difi**3
2136 <         else
2137 <           difi=0.0
2138 <         endif
2139 <        if (energy_dec) then
2140 <         write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
2141 <      &    i,itheta,rad2deg*thetiii,
2142 <      &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
2143 <      &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
2144 <      &    gloc(itheta+nphi-2,icg)
2145 <         endif
2146 <       enddo
2147 <
2148 6081,6082c4660
2149 < C the distributioni.
2150 < ccc        write (iout,*) thetai,thet_pred_mean
2151 ---
2152 > C the distribution.
2153 6112d4689
2154 < C        write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0
2155 6140d4716
2156 < C       write (iout,*) 'termexp',termexp,termm,termpre,i
2157 6181c4757
2158 <       subroutine ebend(etheta,ethetacnstr)
2159 ---
2160 >       subroutine ebend(etheta)
2161 6200d4775
2162 <       include 'COMMON.TORCNSTR'
2163 6208,6211c4783
2164 < c        print *,i,itype(i-1),itype(i),itype(i-2)
2165 <         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
2166 <      &  .or.itype(i).eq.ntyp1) cycle
2167 < C        print *,i,theta(i)
2168 ---
2169 >         if (itype(i-1).eq.ntyp1) cycle
2170 6223,6224c4795
2171 < C        print *,ethetai
2172 <         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
2173 ---
2174 >         if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
2175 6238a4810
2176 >           ityp1=nthetyp+1
2177 6240d4811
2178 <           ityp1=ithetyp((itype(i-2)))
2179 6245c4816
2180 <         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
2181 ---
2182 >         if (i.lt.nres .and. itype(i).ne.ntyp1) then
2183 6260c4831
2184 <           ityp3=ithetyp((itype(i)))
2185 ---
2186 >           ityp3=nthetyp+1
2187 6310d4880
2188 < C       print *,ethetai
2189 6331d4900
2190 < C        print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
2191 6334,6337d4902
2192 < C        print *,"cosph1", (cosph1(k), k=1,nsingle)
2193 < C        print *,"cosph2", (cosph2(k), k=1,nsingle)
2194 < C        print *,"sinph1", (sinph1(k), k=1,nsingle)
2195 < C        print *,"sinph2", (sinph2(k), k=1,nsingle)
2196 6340d4904
2197 < C        print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
2198 6376d4939
2199 < C        print *,ethetai
2200 6385c4948
2201 <         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
2202 ---
2203 >         gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
2204 6387,6415d4949
2205 < C now constrains
2206 <       ethetacnstr=0.0d0
2207 < C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
2208 <       do i=ithetaconstr_start,ithetaconstr_end
2209 <         itheta=itheta_constr(i)
2210 <         thetiii=theta(itheta)
2211 <         difi=pinorm(thetiii-theta_constr0(i))
2212 <         if (difi.gt.theta_drange(i)) then
2213 <           difi=difi-theta_drange(i)
2214 <           ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
2215 <           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
2216 <      &    +for_thet_constr(i)*difi**3
2217 <         else if (difi.lt.-drange(i)) then
2218 <           difi=difi+drange(i)
2219 <           ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
2220 <           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
2221 <      &    +for_thet_constr(i)*difi**3
2222 <         else
2223 <           difi=0.0
2224 <         endif
2225 <        if (energy_dec) then
2226 <         write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
2227 <      &    i,itheta,rad2deg*thetiii,
2228 <      &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
2229 <      &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
2230 <      &    gloc(itheta+nphi-2,icg)
2231 <         endif
2232 <       enddo
2233 <
2234 7180,7182c5714,5716
2235 <      &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
2236 <         itori=itortyp(itype(i-2))
2237 <         itori1=itortyp(itype(i-1))
2238 ---
2239 >      &      .or. itype(i).eq.ntyp1) cycle
2240 >       itori=itortyp(itype(i-2))
2241 >       itori1=itortyp(itype(i-1))
2242 7235,7236c5769,5770
2243 <           edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2244 <           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2245 ---
2246 >           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2247 >           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
2248 7239,7240c5773,5774
2249 <           edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4
2250 <           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2251 ---
2252 >           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2253 >           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
2254 7276,7284c5810,5811
2255 < C ANY TWO ARE DUMMY ATOMS in row CYCLE
2256 < c        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
2257 < c     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
2258 < c     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
2259 <         if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
2260 <      &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
2261 < C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
2262 < C For introducing the NH3+ and COO- group please check the etor_d for reference
2263 < C and guidance
2264 ---
2265 >         if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
2266 >      &       .or. itype(i).eq.ntyp1) cycle
2267 7346,7347c5873,5874
2268 <           edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2269 <           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2270 ---
2271 >           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2272 >           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
2273 7350,7351c5877,5878
2274 <           edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2275 <           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2276 ---
2277 >           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2278 >           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
2279 7355,7360c5882,5884
2280 <        if (energy_dec) then
2281 <         write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc",
2282 <      &    i,itori,rad2deg*phii,
2283 <      &    rad2deg*phi0(i),  rad2deg*drange(i),
2284 <      &    rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg)
2285 <         endif
2286 ---
2287 > cd        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
2288 > cd     &    rad2deg*phi0(i),  rad2deg*drange(i),
2289 > cd     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
2290 7388,7396c5912,5913
2291 < C ANY TWO ARE DUMMY ATOMS in row CYCLE
2292 < C        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
2293 < C     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or.
2294 < C     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))  .or.
2295 < C     &      ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle
2296 <          if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
2297 <      &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
2298 <      &  (itype(i+1).eq.ntyp1)) cycle
2299 < C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
2300 ---
2301 >         if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
2302 >      &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
2303 7406,7414d5922
2304 < C Iblock=2 Proline type
2305 < C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT
2306 < C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO-
2307 < C        if (itype(i+1).eq.ntyp1) iblock=3
2308 < C The problem of NH3+ group can be resolved by adding new parameters please note if there
2309 < C IS or IS NOT need for this
2310 < C IF Yes uncomment below and add to parmread.F appropriate changes and to v1cij and so on
2311 < C        is (itype(i-3).eq.ntyp1) ntblock=2
2312 < C        ntblock is N-terminal blocking group
2313 7418,7420d5925
2314 < C Example of changes for NH3+ blocking group
2315 < C       do j=1,ntermd_1(itori,itori1,itori2,iblock,ntblock)
2316 < C          v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock)
2317 7458,7703d5962
2318 < C----------------------------------------------------------------------------------
2319 < C The rigorous attempt to derive energy function
2320 <       subroutine etor_kcc(etors,edihcnstr)
2321 <       implicit real*8 (a-h,o-z)
2322 <       include 'DIMENSIONS'
2323 <       include 'COMMON.VAR'
2324 <       include 'COMMON.GEO'
2325 <       include 'COMMON.LOCAL'
2326 <       include 'COMMON.TORSION'
2327 <       include 'COMMON.INTERACT'
2328 <       include 'COMMON.DERIV'
2329 <       include 'COMMON.CHAIN'
2330 <       include 'COMMON.NAMES'
2331 <       include 'COMMON.IOUNITS'
2332 <       include 'COMMON.FFIELD'
2333 <       include 'COMMON.TORCNSTR'
2334 <       include 'COMMON.CONTROL'
2335 <       logical lprn
2336 < c      double precision thybt1(maxtermkcc),thybt2(maxtermkcc)
2337 < C Set lprn=.true. for debugging
2338 <       lprn=.false.
2339 < c     lprn=.true.
2340 < C      print *,"wchodze kcc"
2341 <       if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode
2342 <       if (tor_mode.ne.2) then
2343 <       etors=0.0D0
2344 <       endif
2345 <       do i=iphi_start,iphi_end
2346 < C ANY TWO ARE DUMMY ATOMS in row CYCLE
2347 < c        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
2348 < c     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
2349 < c     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
2350 <         if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
2351 <      &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
2352 <         itori=itortyp_kcc(itype(i-2))
2353 <         itori1=itortyp_kcc(itype(i-1))
2354 <         phii=phi(i)
2355 <         glocig=0.0D0
2356 <         glocit1=0.0d0
2357 <         glocit2=0.0d0
2358 <         sumnonchebyshev=0.0d0
2359 <         sumchebyshev=0.0d0
2360 < C to avoid multiple devision by 2
2361 < c        theti22=0.5d0*theta(i)
2362 < C theta 12 is the theta_1 /2
2363 < C theta 22 is theta_2 /2
2364 < c        theti12=0.5d0*theta(i-1)
2365 < C and appropriate sinus function
2366 <         sinthet1=dsin(theta(i-1))
2367 <         sinthet2=dsin(theta(i))
2368 <         costhet1=dcos(theta(i-1))
2369 <         costhet2=dcos(theta(i))
2370 < c Cosines of halves thetas
2371 <         costheti12=0.5d0*(1.0d0+costhet1)
2372 <         costheti22=0.5d0*(1.0d0+costhet2)
2373 < C to speed up lets store its mutliplication
2374 <         sint1t2=sinthet2*sinthet1
2375 <         sint1t2n=1.0d0
2376 < C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma)
2377 < C +d_n*sin(n*gamma)) *
2378 < C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2)))
2379 < C we have two sum 1) Non-Chebyshev which is with n and gamma
2380 <         etori=0.0d0
2381 <         do j=1,nterm_kcc(itori,itori1)
2382 <
2383 <           nval=nterm_kcc_Tb(itori,itori1)
2384 <           v1ij=v1_kcc(j,itori,itori1)
2385 <           v2ij=v2_kcc(j,itori,itori1)
2386 < c          write (iout,*) "i",i," j",j," v1",v1ij," v2",v2ij
2387 < C v1ij is c_n and d_n in euation above
2388 <           cosphi=dcos(j*phii)
2389 <           sinphi=dsin(j*phii)
2390 <           sint1t2n1=sint1t2n
2391 <           sint1t2n=sint1t2n*sint1t2
2392 <           sumth1tyb1=tschebyshev(1,nval,v11_chyb(1,j,itori,itori1),
2393 <      &        costheti12)
2394 <           gradth1tyb1=-0.5d0*sinthet1*gradtschebyshev(0,nval-1,
2395 <      &        v11_chyb(1,j,itori,itori1),costheti12)
2396 < c          write (iout,*) "v11",(v11_chyb(k,j,itori,itori1),k=1,nval),
2397 < c     &      " sumth1tyb1",sumth1tyb1," gradth1tyb1",gradth1tyb1
2398 <           sumth2tyb1=tschebyshev(1,nval,v21_chyb(1,j,itori,itori1),
2399 <      &        costheti22)
2400 <           gradth2tyb1=-0.5d0*sinthet2*gradtschebyshev(0,nval-1,
2401 <      &        v21_chyb(1,j,itori,itori1),costheti22)
2402 < c          write (iout,*) "v21",(v21_chyb(k,j,itori,itori1),k=1,nval),
2403 < c     &      " sumth2tyb1",sumth2tyb1," gradth2tyb1",gradth2tyb1
2404 <           sumth1tyb2=tschebyshev(1,nval,v12_chyb(1,j,itori,itori1),
2405 <      &        costheti12)
2406 <           gradth1tyb2=-0.5d0*sinthet1*gradtschebyshev(0,nval-1,
2407 <      &        v12_chyb(1,j,itori,itori1),costheti12)
2408 < c          write (iout,*) "v12",(v12_chyb(k,j,itori,itori1),k=1,nval),
2409 < c     &      " sumth1tyb2",sumth1tyb2," gradth1tyb2",gradth1tyb2
2410 <           sumth2tyb2=tschebyshev(1,nval,v22_chyb(1,j,itori,itori1),
2411 <      &        costheti22)
2412 <           gradth2tyb2=-0.5d0*sinthet2*gradtschebyshev(0,nval-1,
2413 <      &        v22_chyb(1,j,itori,itori1),costheti22)
2414 < c          write (iout,*) "v22",(v22_chyb(k,j,itori,itori1),k=1,nval),
2415 < c     &      " sumth2tyb2",sumth2tyb2," gradth2tyb2",gradth2tyb2
2416 < C          etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi)
2417 < C          if (energy_dec) etors_ii=etors_ii+
2418 < C     &                v1ij*cosphi+v2ij*sinphi
2419 < C glocig is the gradient local i site in gamma
2420 <           actval1=v1ij*cosphi*(1.0d0+sumth1tyb1+sumth2tyb1)
2421 <           actval2=v2ij*sinphi*(1.0d0+sumth1tyb2+sumth2tyb2)
2422 <           etori=etori+sint1t2n*(actval1+actval2)
2423 <           glocig=glocig+
2424 <      &        j*sint1t2n*(v2ij*cosphi*(1.0d0+sumth1tyb2+sumth2tyb2)
2425 <      &        -v1ij*sinphi*(1.0d0+sumth1tyb1+sumth2tyb1))
2426 < C now gradient over theta_1
2427 <           glocit1=glocit1+
2428 <      &       j*sint1t2n1*costhet1*sinthet2*(actval1+actval2)+
2429 <      &       sint1t2n*(v1ij*cosphi*gradth1tyb1+v2ij*sinphi*gradth1tyb2)
2430 <           glocit2=glocit2+
2431 <      &       j*sint1t2n1*sinthet1*costhet2*(actval1+actval2)+
2432 <      &       sint1t2n*(v1ij*cosphi*gradth2tyb1+v2ij*sinphi*gradth2tyb2)
2433 <
2434 < C now the Czebyshev polinominal sum
2435 < c        do k=1,nterm_kcc_Tb(itori,itori1)
2436 < c         thybt1(k)=v1_chyb(k,j,itori,itori1)
2437 < c         thybt2(k)=v2_chyb(k,j,itori,itori1)
2438 < C         thybt1(k)=0.0
2439 < C         thybt2(k)=0.0
2440 < c        enddo
2441 < C        print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2,
2442 < C     &         gradtschebyshev
2443 < C     &         (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1),
2444 < C     &         dcos(theti22)**2),
2445 < C     &         dsin(theti22)
2446 <
2447 < C now overal sumation
2448 < C         print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb
2449 <         enddo ! j
2450 <         etors=etors+etori
2451 < C derivative over gamma
2452 <         gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
2453 < C derivative over theta1
2454 <         gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1
2455 < C now derivative over theta2
2456 <         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2
2457 <         if (lprn)
2458 <      &    write (iout,*) i-2,i-1,itype(i-2),itype(i-1),itori,itori1,
2459 <      &       theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori
2460 <       enddo
2461 < C        gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
2462 < ! 6/20/98 - dihedral angle constraints
2463 <       if (tor_mode.ne.2) then
2464 <       edihcnstr=0.0d0
2465 < c      do i=1,ndih_constr
2466 <       do i=idihconstr_start,idihconstr_end
2467 <         itori=idih_constr(i)
2468 <         phii=phi(itori)
2469 <         difi=pinorm(phii-phi0(i))
2470 <         if (difi.gt.drange(i)) then
2471 <           difi=difi-drange(i)
2472 <           edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2473 <           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2474 <         else if (difi.lt.-drange(i)) then
2475 <           difi=difi+drange(i)
2476 <           edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2477 <           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2478 <         else
2479 <           difi=0.0
2480 <         endif
2481 <        enddo
2482 <        endif
2483 <       return
2484 <       end
2485 <
2486 < C The rigorous attempt to derive energy function
2487 <       subroutine ebend_kcc(etheta,ethetacnstr)
2488 <
2489 <       implicit real*8 (a-h,o-z)
2490 <       include 'DIMENSIONS'
2491 <       include 'COMMON.VAR'
2492 <       include 'COMMON.GEO'
2493 <       include 'COMMON.LOCAL'
2494 <       include 'COMMON.TORSION'
2495 <       include 'COMMON.INTERACT'
2496 <       include 'COMMON.DERIV'
2497 <       include 'COMMON.CHAIN'
2498 <       include 'COMMON.NAMES'
2499 <       include 'COMMON.IOUNITS'
2500 <       include 'COMMON.FFIELD'
2501 <       include 'COMMON.TORCNSTR'
2502 <       include 'COMMON.CONTROL'
2503 <       logical lprn
2504 <       double precision thybt1(maxtermkcc)
2505 < C Set lprn=.true. for debugging
2506 <       lprn=.false.
2507 < c     lprn=.true.
2508 < C      print *,"wchodze kcc"
2509 <       if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode
2510 <       if (tor_mode.ne.2) etheta=0.0D0
2511 <       do i=ithet_start,ithet_end
2512 < c        print *,i,itype(i-1),itype(i),itype(i-2)
2513 <         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
2514 <      &  .or.itype(i).eq.ntyp1) cycle
2515 <          iti=itortyp_kcc(itype(i-1))
2516 <         sinthet=dsin(theta(i)/2.0d0)
2517 <         costhet=dcos(theta(i)/2.0d0)
2518 <          do j=1,nbend_kcc_Tb(iti)
2519 <           thybt1(j)=v1bend_chyb(j,iti)
2520 <          enddo
2521 <          sumth1thyb=tschebyshev
2522 <      &         (1,nbend_kcc_Tb(iti),thybt1(1),costhet)
2523 <         if (lprn) write (iout,*) i-1,itype(i-1),iti,theta(i)*rad2deg,
2524 <      &    sumth1thyb
2525 <         ihelp=nbend_kcc_Tb(iti)-1
2526 <         gradthybt1=gradtschebyshev
2527 <      &         (0,ihelp,thybt1(1),costhet)
2528 <         etheta=etheta+sumth1thyb
2529 < C        print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
2530 <         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*
2531 <      &   gradthybt1*sinthet*(-0.5d0)
2532 <       enddo
2533 <       if (tor_mode.ne.2) then
2534 <       ethetacnstr=0.0d0
2535 < C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
2536 <       do i=ithetaconstr_start,ithetaconstr_end
2537 <         itheta=itheta_constr(i)
2538 <         thetiii=theta(itheta)
2539 <         difi=pinorm(thetiii-theta_constr0(i))
2540 <         if (difi.gt.theta_drange(i)) then
2541 <           difi=difi-theta_drange(i)
2542 <           ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
2543 <           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
2544 <      &    +for_thet_constr(i)*difi**3
2545 <         else if (difi.lt.-drange(i)) then
2546 <           difi=difi+drange(i)
2547 <           ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
2548 <           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
2549 <      &    +for_thet_constr(i)*difi**3
2550 <         else
2551 <           difi=0.0
2552 <         endif
2553 <        if (energy_dec) then
2554 <         write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
2555 <      &    i,itheta,rad2deg*thetiii,
2556 <      &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
2557 <      &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
2558 <      &    gloc(itheta+nphi-2,icg)
2559 <         endif
2560 <       enddo
2561 <       endif
2562 <       return
2563 <       end
2564 7845d6103
2565 <       include 'COMMON.SHIELD'
2566 8264d6521
2567 <       include 'COMMON.SHIELD'
2568 8567d6823
2569 < CC     &            *fac_shield(i)**2*fac_shield(j)**2
2570 8687,8688d6942
2571 <       include 'COMMON.SHIELD'
2572 <       include 'COMMON.CONTROL'
2573 8692d6945
2574 < C      print *,"wchodze",fac_shield(i),shield_mode
2575 8701,8702d6953
2576 < C*
2577 < C     & fac_shield(i)**2*fac_shield(j)**2
2578 8766d7016
2579 < C      print *,ekont,ees,i,k
2580 8768,8848d7017
2581 < C now gradient over shielding
2582 < C      return
2583 <       if (shield_mode.gt.0) then
2584 <        j=ees0plist(jj,i)
2585 <        l=ees0plist(kk,k)
2586 < C        print *,i,j,fac_shield(i),fac_shield(j),
2587 < C     &fac_shield(k),fac_shield(l)
2588 <         if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
2589 <      &      (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then
2590 <           do ilist=1,ishield_list(i)
2591 <            iresshield=shield_list(ilist,i)
2592 <            do m=1,3
2593 <            rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i)
2594 < C     &      *2.0
2595 <            gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2596 <      &              rlocshield
2597 <      & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i)
2598 <             gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2599 <      &+rlocshield
2600 <            enddo
2601 <           enddo
2602 <           do ilist=1,ishield_list(j)
2603 <            iresshield=shield_list(ilist,j)
2604 <            do m=1,3
2605 <            rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j)
2606 < C     &     *2.0
2607 <            gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2608 <      &              rlocshield
2609 <      & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j)
2610 <            gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2611 <      &     +rlocshield
2612 <            enddo
2613 <           enddo
2614 <
2615 <           do ilist=1,ishield_list(k)
2616 <            iresshield=shield_list(ilist,k)
2617 <            do m=1,3
2618 <            rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k)
2619 < C     &     *2.0
2620 <            gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2621 <      &              rlocshield
2622 <      & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k)
2623 <            gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2624 <      &     +rlocshield
2625 <            enddo
2626 <           enddo
2627 <           do ilist=1,ishield_list(l)
2628 <            iresshield=shield_list(ilist,l)
2629 <            do m=1,3
2630 <            rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l)
2631 < C     &     *2.0
2632 <            gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2633 <      &              rlocshield
2634 <      & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l)
2635 <            gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2636 <      &     +rlocshield
2637 <            enddo
2638 <           enddo
2639 < C          print *,gshieldx(m,iresshield)
2640 <           do m=1,3
2641 <             gshieldc_ec(m,i)=gshieldc_ec(m,i)+
2642 <      &              grad_shield(m,i)*ehbcorr/fac_shield(i)
2643 <             gshieldc_ec(m,j)=gshieldc_ec(m,j)+
2644 <      &              grad_shield(m,j)*ehbcorr/fac_shield(j)
2645 <             gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+
2646 <      &              grad_shield(m,i)*ehbcorr/fac_shield(i)
2647 <             gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+
2648 <      &              grad_shield(m,j)*ehbcorr/fac_shield(j)
2649 <
2650 <             gshieldc_ec(m,k)=gshieldc_ec(m,k)+
2651 <      &              grad_shield(m,k)*ehbcorr/fac_shield(k)
2652 <             gshieldc_ec(m,l)=gshieldc_ec(m,l)+
2653 <      &              grad_shield(m,l)*ehbcorr/fac_shield(l)
2654 <             gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+
2655 <      &              grad_shield(m,k)*ehbcorr/fac_shield(k)
2656 <             gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+
2657 <      &              grad_shield(m,l)*ehbcorr/fac_shield(l)
2658 <
2659 <            enddo
2660 <       endif
2661 <       endif
2662 8867c7036
2663 <       iti1 = itortyp(itype(i+1))
2664 ---
2665 >       iti1 = itype2loc(itype(i+1))
2666 9531c7700
2667 <       call transpose2(EE(1,1,k),auxmat(1,1))
2668 ---
2669 >       call transpose2(EE(1,1,itk),auxmat(1,1))
2670 9612c7781
2671 <         call transpose2(EE(1,1,l),auxmat(1,1))
2672 ---
2673 >         call transpose2(EE(1,1,itl),auxmat(1,1))
2674 9685c7854
2675 <         call transpose2(EE(1,1,j),auxmat(1,1))
2676 ---
2677 >         call transpose2(EE(1,1,itj),auxmat(1,1))
2678 9773c7942
2679 <         gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2)
2680 ---
2681 >         gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2)
2682 9775c7944
2683 <         gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2)
2684 ---
2685 >         gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2)
2686 10080c8249
2687 < C       o| o |         | o |o                                                  C
2688 ---
2689 > C       o| o |         | o |o                                                  C
2690 10084,10085c8253,8254
2691 < C       i             i                                                        C
2692 < C                                                                              C
2693 ---
2694 > C       i             i                                                        C
2695 > C                                                                              C
2696 10256c8425
2697 < C                                                                              C
2698 ---
2699 > C                                                                              C
2700 10259c8428
2701 < C          o             o                                                     C
2702 ---
2703 > C          o             o                                                     C
2704 10272c8441
2705 <       iti=itortyp(itype(i))
2706 ---
2707 >       iti=itype2loc(itype(i))
2708 10292c8461
2709 <       call transpose2(EE(1,1,k),auxmat(1,1))
2710 ---
2711 >       call transpose2(EE(1,1,itk),auxmat(1,1))
2712 10373c8542
2713 < C                                                                              C
2714 ---
2715 > C                                                                              C
2716 10381c8550
2717 < C      \ /   \       \ /   \                                                   C
2718 ---
2719 > C      \ /   \       \ /   \                                                   C
2720 10384c8553
2721 < C                                                                              C
2722 ---
2723 > C                                                                              C
2724 11082,11656d9250
2725 <       return
2726 <       end
2727 < CCC----------------------------------------------
2728 <       subroutine Eliptransfer(eliptran)
2729 <       implicit real*8 (a-h,o-z)
2730 <       include 'DIMENSIONS'
2731 <       include 'COMMON.GEO'
2732 <       include 'COMMON.VAR'
2733 <       include 'COMMON.LOCAL'
2734 <       include 'COMMON.CHAIN'
2735 <       include 'COMMON.DERIV'
2736 <       include 'COMMON.NAMES'
2737 <       include 'COMMON.INTERACT'
2738 <       include 'COMMON.IOUNITS'
2739 <       include 'COMMON.CALC'
2740 <       include 'COMMON.CONTROL'
2741 <       include 'COMMON.SPLITELE'
2742 <       include 'COMMON.SBRIDGE'
2743 < C this is done by Adasko
2744 < C      print *,"wchodze"
2745 < C structure of box:
2746 < C      water
2747 < C--bordliptop-- buffore starts
2748 < C--bufliptop--- here true lipid starts
2749 < C      lipid
2750 < C--buflipbot--- lipid ends buffore starts
2751 < C--bordlipbot--buffore ends
2752 <       eliptran=0.0
2753 <       do i=ilip_start,ilip_end
2754 < C       do i=1,1
2755 <         if (itype(i).eq.ntyp1) cycle
2756 <
2757 <         positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
2758 <         if (positi.le.0) positi=positi+boxzsize
2759 < C        print *,i
2760 < C first for peptide groups
2761 < c for each residue check if it is in lipid or lipid water border area
2762 <        if ((positi.gt.bordlipbot)
2763 <      &.and.(positi.lt.bordliptop)) then
2764 < C the energy transfer exist
2765 <         if (positi.lt.buflipbot) then
2766 < C what fraction I am in
2767 <          fracinbuf=1.0d0-
2768 <      &        ((positi-bordlipbot)/lipbufthick)
2769 < C lipbufthick is thickenes of lipid buffore
2770 <          sslip=sscalelip(fracinbuf)
2771 <          ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
2772 <          eliptran=eliptran+sslip*pepliptran
2773 <          gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
2774 <          gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
2775 < C         gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
2776 <
2777 < C        print *,"doing sccale for lower part"
2778 < C         print *,i,sslip,fracinbuf,ssgradlip
2779 <         elseif (positi.gt.bufliptop) then
2780 <          fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
2781 <          sslip=sscalelip(fracinbuf)
2782 <          ssgradlip=sscagradlip(fracinbuf)/lipbufthick
2783 <          eliptran=eliptran+sslip*pepliptran
2784 <          gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
2785 <          gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
2786 < C         gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
2787 < C          print *, "doing sscalefor top part"
2788 < C         print *,i,sslip,fracinbuf,ssgradlip
2789 <         else
2790 <          eliptran=eliptran+pepliptran
2791 < C         print *,"I am in true lipid"
2792 <         endif
2793 < C       else
2794 < C       eliptran=elpitran+0.0 ! I am in water
2795 <        endif
2796 <        enddo
2797 < C       print *, "nic nie bylo w lipidzie?"
2798 < C now multiply all by the peptide group transfer factor
2799 < C       eliptran=eliptran*pepliptran
2800 < C now the same for side chains
2801 < CV       do i=1,1
2802 <        do i=ilip_start,ilip_end
2803 <         if (itype(i).eq.ntyp1) cycle
2804 <         positi=(mod(c(3,i+nres),boxzsize))
2805 <         if (positi.le.0) positi=positi+boxzsize
2806 < C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
2807 < c for each residue check if it is in lipid or lipid water border area
2808 < C       respos=mod(c(3,i+nres),boxzsize)
2809 < C       print *,positi,bordlipbot,buflipbot
2810 <        if ((positi.gt.bordlipbot)
2811 <      & .and.(positi.lt.bordliptop)) then
2812 < C the energy transfer exist
2813 <         if (positi.lt.buflipbot) then
2814 <          fracinbuf=1.0d0-
2815 <      &     ((positi-bordlipbot)/lipbufthick)
2816 < C lipbufthick is thickenes of lipid buffore
2817 <          sslip=sscalelip(fracinbuf)
2818 <          ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
2819 <          eliptran=eliptran+sslip*liptranene(itype(i))
2820 <          gliptranx(3,i)=gliptranx(3,i)
2821 <      &+ssgradlip*liptranene(itype(i))
2822 <          gliptranc(3,i-1)= gliptranc(3,i-1)
2823 <      &+ssgradlip*liptranene(itype(i))
2824 < C         print *,"doing sccale for lower part"
2825 <         elseif (positi.gt.bufliptop) then
2826 <          fracinbuf=1.0d0-
2827 <      &((bordliptop-positi)/lipbufthick)
2828 <          sslip=sscalelip(fracinbuf)
2829 <          ssgradlip=sscagradlip(fracinbuf)/lipbufthick
2830 <          eliptran=eliptran+sslip*liptranene(itype(i))
2831 <          gliptranx(3,i)=gliptranx(3,i)
2832 <      &+ssgradlip*liptranene(itype(i))
2833 <          gliptranc(3,i-1)= gliptranc(3,i-1)
2834 <      &+ssgradlip*liptranene(itype(i))
2835 < C          print *, "doing sscalefor top part",sslip,fracinbuf
2836 <         else
2837 <          eliptran=eliptran+liptranene(itype(i))
2838 < C         print *,"I am in true lipid"
2839 <         endif
2840 <         endif ! if in lipid or buffor
2841 < C       else
2842 < C       eliptran=elpitran+0.0 ! I am in water
2843 <        enddo
2844 <        return
2845 <        end
2846 < C---------------------------------------------------------
2847 < C AFM soubroutine for constant force
2848 <        subroutine AFMforce(Eafmforce)
2849 <        implicit real*8 (a-h,o-z)
2850 <       include 'DIMENSIONS'
2851 <       include 'COMMON.GEO'
2852 <       include 'COMMON.VAR'
2853 <       include 'COMMON.LOCAL'
2854 <       include 'COMMON.CHAIN'
2855 <       include 'COMMON.DERIV'
2856 <       include 'COMMON.NAMES'
2857 <       include 'COMMON.INTERACT'
2858 <       include 'COMMON.IOUNITS'
2859 <       include 'COMMON.CALC'
2860 <       include 'COMMON.CONTROL'
2861 <       include 'COMMON.SPLITELE'
2862 <       include 'COMMON.SBRIDGE'
2863 <       real*8 diffafm(3)
2864 <       dist=0.0d0
2865 <       Eafmforce=0.0d0
2866 <       do i=1,3
2867 <       diffafm(i)=c(i,afmend)-c(i,afmbeg)
2868 <       dist=dist+diffafm(i)**2
2869 <       enddo
2870 <       dist=dsqrt(dist)
2871 <       Eafmforce=-forceAFMconst*(dist-distafminit)
2872 <       do i=1,3
2873 <       gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist
2874 <       gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist
2875 <       enddo
2876 < C      print *,'AFM',Eafmforce
2877 <       return
2878 <       end
2879 < C---------------------------------------------------------
2880 < C AFM subroutine with pseudoconstant velocity
2881 <        subroutine AFMvel(Eafmforce)
2882 <        implicit real*8 (a-h,o-z)
2883 <       include 'DIMENSIONS'
2884 <       include 'COMMON.GEO'
2885 <       include 'COMMON.VAR'
2886 <       include 'COMMON.LOCAL'
2887 <       include 'COMMON.CHAIN'
2888 <       include 'COMMON.DERIV'
2889 <       include 'COMMON.NAMES'
2890 <       include 'COMMON.INTERACT'
2891 <       include 'COMMON.IOUNITS'
2892 <       include 'COMMON.CALC'
2893 <       include 'COMMON.CONTROL'
2894 <       include 'COMMON.SPLITELE'
2895 <       include 'COMMON.SBRIDGE'
2896 <       real*8 diffafm(3)
2897 < C Only for check grad COMMENT if not used for checkgrad
2898 < C      totT=3.0d0
2899 < C--------------------------------------------------------
2900 < C      print *,"wchodze"
2901 <       dist=0.0d0
2902 <       Eafmforce=0.0d0
2903 <       do i=1,3
2904 <       diffafm(i)=c(i,afmend)-c(i,afmbeg)
2905 <       dist=dist+diffafm(i)**2
2906 <       enddo
2907 <       dist=dsqrt(dist)
2908 <       Eafmforce=0.5d0*forceAFMconst
2909 <      & *(distafminit+totTafm*velAFMconst-dist)**2
2910 < C      Eafmforce=-forceAFMconst*(dist-distafminit)
2911 <       do i=1,3
2912 <       gradafm(i,afmend-1)=-forceAFMconst*
2913 <      &(distafminit+totTafm*velAFMconst-dist)
2914 <      &*diffafm(i)/dist
2915 <       gradafm(i,afmbeg-1)=forceAFMconst*
2916 <      &(distafminit+totTafm*velAFMconst-dist)
2917 <      &*diffafm(i)/dist
2918 <       enddo
2919 < C      print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
2920 <       return
2921 <       end
2922 < C-----------------------------------------------------------
2923 < C first for shielding is setting of function of side-chains
2924 <        subroutine set_shield_fac
2925 <       implicit real*8 (a-h,o-z)
2926 <       include 'DIMENSIONS'
2927 <       include 'COMMON.CHAIN'
2928 <       include 'COMMON.DERIV'
2929 <       include 'COMMON.IOUNITS'
2930 <       include 'COMMON.SHIELD'
2931 <       include 'COMMON.INTERACT'
2932 < C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
2933 <       double precision div77_81/0.974996043d0/,
2934 <      &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
2935 <
2936 < C the vector between center of side_chain and peptide group
2937 <        double precision pep_side(3),long,side_calf(3),
2938 <      &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
2939 <      &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
2940 < C the line belowe needs to be changed for FGPROC>1
2941 <       do i=1,nres-1
2942 <       if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
2943 <       ishield_list(i)=0
2944 < Cif there two consequtive dummy atoms there is no peptide group between them
2945 < C the line below has to be changed for FGPROC>1
2946 <       VolumeTotal=0.0
2947 <       do k=1,nres
2948 <        if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
2949 <        dist_pep_side=0.0
2950 <        dist_side_calf=0.0
2951 <        do j=1,3
2952 < C first lets set vector conecting the ithe side-chain with kth side-chain
2953 <       pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
2954 < C      pep_side(j)=2.0d0
2955 < C and vector conecting the side-chain with its proper calfa
2956 <       side_calf(j)=c(j,k+nres)-c(j,k)
2957 < C      side_calf(j)=2.0d0
2958 <       pept_group(j)=c(j,i)-c(j,i+1)
2959 < C lets have their lenght
2960 <       dist_pep_side=pep_side(j)**2+dist_pep_side
2961 <       dist_side_calf=dist_side_calf+side_calf(j)**2
2962 <       dist_pept_group=dist_pept_group+pept_group(j)**2
2963 <       enddo
2964 <        dist_pep_side=dsqrt(dist_pep_side)
2965 <        dist_pept_group=dsqrt(dist_pept_group)
2966 <        dist_side_calf=dsqrt(dist_side_calf)
2967 <       do j=1,3
2968 <         pep_side_norm(j)=pep_side(j)/dist_pep_side
2969 <         side_calf_norm(j)=dist_side_calf
2970 <       enddo
2971 < C now sscale fraction
2972 <        sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
2973 < C       print *,buff_shield,"buff"
2974 < C now sscale
2975 <         if (sh_frac_dist.le.0.0) cycle
2976 < C If we reach here it means that this side chain reaches the shielding sphere
2977 < C Lets add him to the list for gradient
2978 <         ishield_list(i)=ishield_list(i)+1
2979 < C ishield_list is a list of non 0 side-chain that contribute to factor gradient
2980 < C this list is essential otherwise problem would be O3
2981 <         shield_list(ishield_list(i),i)=k
2982 < C Lets have the sscale value
2983 <         if (sh_frac_dist.gt.1.0) then
2984 <          scale_fac_dist=1.0d0
2985 <          do j=1,3
2986 <          sh_frac_dist_grad(j)=0.0d0
2987 <          enddo
2988 <         else
2989 <          scale_fac_dist=-sh_frac_dist*sh_frac_dist
2990 <      &                   *(2.0*sh_frac_dist-3.0d0)
2991 <          fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
2992 <      &                  /dist_pep_side/buff_shield*0.5
2993 < C remember for the final gradient multiply sh_frac_dist_grad(j)
2994 < C for side_chain by factor -2 !
2995 <          do j=1,3
2996 <          sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
2997 < C         print *,"jestem",scale_fac_dist,fac_help_scale,
2998 < C     &                    sh_frac_dist_grad(j)
2999 <          enddo
3000 <         endif
3001 < C        if ((i.eq.3).and.(k.eq.2)) then
3002 < C        print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
3003 < C     & ,"TU"
3004 < C        endif
3005 <
3006 < C this is what is now we have the distance scaling now volume...
3007 <       short=short_r_sidechain(itype(k))
3008 <       long=long_r_sidechain(itype(k))
3009 <       costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
3010 < C now costhet_grad
3011 < C       costhet=0.0d0
3012 <        costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
3013 < C       costhet_fac=0.0d0
3014 <        do j=1,3
3015 <          costhet_grad(j)=costhet_fac*pep_side(j)
3016 <        enddo
3017 < C remember for the final gradient multiply costhet_grad(j)
3018 < C for side_chain by factor -2 !
3019 < C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
3020 < C pep_side0pept_group is vector multiplication
3021 <       pep_side0pept_group=0.0
3022 <       do j=1,3
3023 <       pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
3024 <       enddo
3025 <       cosalfa=(pep_side0pept_group/
3026 <      & (dist_pep_side*dist_side_calf))
3027 <       fac_alfa_sin=1.0-cosalfa**2
3028 <       fac_alfa_sin=dsqrt(fac_alfa_sin)
3029 <       rkprim=fac_alfa_sin*(long-short)+short
3030 < C now costhet_grad
3031 <        cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
3032 <        cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
3033 <
3034 <        do j=1,3
3035 <          cosphi_grad_long(j)=cosphi_fac*pep_side(j)
3036 <      &+cosphi**3*0.5/dist_pep_side**2*(-rkprim)
3037 <      &*(long-short)/fac_alfa_sin*cosalfa/
3038 <      &((dist_pep_side*dist_side_calf))*
3039 <      &((side_calf(j))-cosalfa*
3040 <      &((pep_side(j)/dist_pep_side)*dist_side_calf))
3041 <
3042 <         cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim)
3043 <      &*(long-short)/fac_alfa_sin*cosalfa
3044 <      &/((dist_pep_side*dist_side_calf))*
3045 <      &(pep_side(j)-
3046 <      &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
3047 <        enddo
3048 <
3049 <       VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
3050 <      &                    /VSolvSphere_div
3051 <      &                    *wshield
3052 < C now the gradient...
3053 < C grad_shield is gradient of Calfa for peptide groups
3054 < C      write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
3055 < C     &               costhet,cosphi
3056 < C       write(iout,*) "cosphi_compon",i,k,pep_side0pept_group,
3057 < C     & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k)
3058 <       do j=1,3
3059 <       grad_shield(j,i)=grad_shield(j,i)
3060 < C gradient po skalowaniu
3061 <      &                +(sh_frac_dist_grad(j)
3062 < C  gradient po costhet
3063 <      &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
3064 <      &-scale_fac_dist*(cosphi_grad_long(j))
3065 <      &/(1.0-cosphi) )*div77_81
3066 <      &*VofOverlap
3067 < C grad_shield_side is Cbeta sidechain gradient
3068 <       grad_shield_side(j,ishield_list(i),i)=
3069 <      &        (sh_frac_dist_grad(j)*-2.0d0
3070 <      &       +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
3071 <      &       +scale_fac_dist*(cosphi_grad_long(j))
3072 <      &        *2.0d0/(1.0-cosphi))
3073 <      &        *div77_81*VofOverlap
3074 <
3075 <        grad_shield_loc(j,ishield_list(i),i)=
3076 <      &   scale_fac_dist*cosphi_grad_loc(j)
3077 <      &        *2.0d0/(1.0-cosphi)
3078 <      &        *div77_81*VofOverlap
3079 <       enddo
3080 <       VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
3081 <       enddo
3082 <       fac_shield(i)=VolumeTotal*div77_81+div4_81
3083 < C      write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
3084 <       enddo
3085 <       return
3086 <       end
3087 < C--------------------------------------------------------------------------
3088 <       double precision function tschebyshev(m,n,x,y)
3089 <       implicit none
3090 <       include "DIMENSIONS"
3091 <       integer i,m,n
3092 <       double precision x(n),y,yy(0:maxvar),aux
3093 < c Tschebyshev polynomial. Note that the first term is omitted
3094 < c m=0: the constant term is included
3095 < c m=1: the constant term is not included
3096 <       yy(0)=1.0d0
3097 <       yy(1)=y
3098 <       do i=2,n
3099 <         yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
3100 <       enddo
3101 <       aux=0.0d0
3102 <       do i=m,n
3103 <         aux=aux+x(i)*yy(i)
3104 <       enddo
3105 <       tschebyshev=aux
3106 <       return
3107 <       end
3108 < C--------------------------------------------------------------------------
3109 <       double precision function gradtschebyshev(m,n,x,y)
3110 <       implicit none
3111 <       include "DIMENSIONS"
3112 <       integer i,m,n
3113 <       double precision x(n+1),y,yy(0:maxvar),aux
3114 < c Tschebyshev polynomial. Note that the first term is omitted
3115 < c m=0: the constant term is included
3116 < c m=1: the constant term is not included
3117 <       yy(0)=1.0d0
3118 <       yy(1)=2.0d0*y
3119 <       do i=2,n
3120 <         yy(i)=2*y*yy(i-1)-yy(i-2)
3121 <       enddo
3122 <       aux=0.0d0
3123 <       do i=m,n
3124 <         aux=aux+x(i+1)*yy(i)*(i+1)
3125 < C        print *, x(i+1),yy(i),i
3126 <       enddo
3127 <       gradtschebyshev=aux
3128 <       return
3129 <       end
3130 < C------------------------------------------------------------------------
3131 < C first for shielding is setting of function of side-chains
3132 <        subroutine set_shield_fac2
3133 <       implicit real*8 (a-h,o-z)
3134 <       include 'DIMENSIONS'
3135 <       include 'COMMON.CHAIN'
3136 <       include 'COMMON.DERIV'
3137 <       include 'COMMON.IOUNITS'
3138 <       include 'COMMON.SHIELD'
3139 <       include 'COMMON.INTERACT'
3140 < C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
3141 <       double precision div77_81/0.974996043d0/,
3142 <      &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
3143 <
3144 < C the vector between center of side_chain and peptide group
3145 <        double precision pep_side(3),long,side_calf(3),
3146 <      &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
3147 <      &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
3148 < C the line belowe needs to be changed for FGPROC>1
3149 <       do i=1,nres-1
3150 <       if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
3151 <       ishield_list(i)=0
3152 < Cif there two consequtive dummy atoms there is no peptide group between them
3153 < C the line below has to be changed for FGPROC>1
3154 <       VolumeTotal=0.0
3155 <       do k=1,nres
3156 <        if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
3157 <        dist_pep_side=0.0
3158 <        dist_side_calf=0.0
3159 <        do j=1,3
3160 < C first lets set vector conecting the ithe side-chain with kth side-chain
3161 <       pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
3162 < C      pep_side(j)=2.0d0
3163 < C and vector conecting the side-chain with its proper calfa
3164 <       side_calf(j)=c(j,k+nres)-c(j,k)
3165 < C      side_calf(j)=2.0d0
3166 <       pept_group(j)=c(j,i)-c(j,i+1)
3167 < C lets have their lenght
3168 <       dist_pep_side=pep_side(j)**2+dist_pep_side
3169 <       dist_side_calf=dist_side_calf+side_calf(j)**2
3170 <       dist_pept_group=dist_pept_group+pept_group(j)**2
3171 <       enddo
3172 <        dist_pep_side=dsqrt(dist_pep_side)
3173 <        dist_pept_group=dsqrt(dist_pept_group)
3174 <        dist_side_calf=dsqrt(dist_side_calf)
3175 <       do j=1,3
3176 <         pep_side_norm(j)=pep_side(j)/dist_pep_side
3177 <         side_calf_norm(j)=dist_side_calf
3178 <       enddo
3179 < C now sscale fraction
3180 <        sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
3181 < C       print *,buff_shield,"buff"
3182 < C now sscale
3183 <         if (sh_frac_dist.le.0.0) cycle
3184 < C If we reach here it means that this side chain reaches the shielding sphere
3185 < C Lets add him to the list for gradient
3186 <         ishield_list(i)=ishield_list(i)+1
3187 < C ishield_list is a list of non 0 side-chain that contribute to factor gradient
3188 < C this list is essential otherwise problem would be O3
3189 <         shield_list(ishield_list(i),i)=k
3190 < C Lets have the sscale value
3191 <         if (sh_frac_dist.gt.1.0) then
3192 <          scale_fac_dist=1.0d0
3193 <          do j=1,3
3194 <          sh_frac_dist_grad(j)=0.0d0
3195 <          enddo
3196 <         else
3197 <          scale_fac_dist=-sh_frac_dist*sh_frac_dist
3198 <      &                   *(2.0d0*sh_frac_dist-3.0d0)
3199 <          fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2)
3200 <      &                  /dist_pep_side/buff_shield*0.5d0
3201 < C remember for the final gradient multiply sh_frac_dist_grad(j)
3202 < C for side_chain by factor -2 !
3203 <          do j=1,3
3204 <          sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
3205 < C         sh_frac_dist_grad(j)=0.0d0
3206 < C         scale_fac_dist=1.0d0
3207 < C         print *,"jestem",scale_fac_dist,fac_help_scale,
3208 < C     &                    sh_frac_dist_grad(j)
3209 <          enddo
3210 <         endif
3211 < C this is what is now we have the distance scaling now volume...
3212 <       short=short_r_sidechain(itype(k))
3213 <       long=long_r_sidechain(itype(k))
3214 <       costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
3215 <       sinthet=short/dist_pep_side*costhet
3216 < C now costhet_grad
3217 < C       costhet=0.6d0
3218 < C       sinthet=0.8
3219 <        costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
3220 < C       sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
3221 < C     &             -short/dist_pep_side**2/costhet)
3222 < C       costhet_fac=0.0d0
3223 <        do j=1,3
3224 <          costhet_grad(j)=costhet_fac*pep_side(j)
3225 <        enddo
3226 < C remember for the final gradient multiply costhet_grad(j)
3227 < C for side_chain by factor -2 !
3228 < C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
3229 < C pep_side0pept_group is vector multiplication
3230 <       pep_side0pept_group=0.0d0
3231 <       do j=1,3
3232 <       pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
3233 <       enddo
3234 <       cosalfa=(pep_side0pept_group/
3235 <      & (dist_pep_side*dist_side_calf))
3236 <       fac_alfa_sin=1.0d0-cosalfa**2
3237 <       fac_alfa_sin=dsqrt(fac_alfa_sin)
3238 <       rkprim=fac_alfa_sin*(long-short)+short
3239 < C      rkprim=short
3240 <
3241 < C now costhet_grad
3242 <        cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
3243 < C       cosphi=0.6
3244 <        cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
3245 <        sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/
3246 <      &      dist_pep_side**2)
3247 < C       sinphi=0.8
3248 <        do j=1,3
3249 <          cosphi_grad_long(j)=cosphi_fac*pep_side(j)
3250 <      &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
3251 <      &*(long-short)/fac_alfa_sin*cosalfa/
3252 <      &((dist_pep_side*dist_side_calf))*
3253 <      &((side_calf(j))-cosalfa*
3254 <      &((pep_side(j)/dist_pep_side)*dist_side_calf))
3255 < C       cosphi_grad_long(j)=0.0d0
3256 <         cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
3257 <      &*(long-short)/fac_alfa_sin*cosalfa
3258 <      &/((dist_pep_side*dist_side_calf))*
3259 <      &(pep_side(j)-
3260 <      &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
3261 < C       cosphi_grad_loc(j)=0.0d0
3262 <        enddo
3263 < C      print *,sinphi,sinthet
3264 <       VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
3265 <      &                    /VSolvSphere_div
3266 < C     &                    *wshield
3267 < C now the gradient...
3268 <       do j=1,3
3269 <       grad_shield(j,i)=grad_shield(j,i)
3270 < C gradient po skalowaniu
3271 <      &                +(sh_frac_dist_grad(j)*VofOverlap
3272 < C  gradient po costhet
3273 <      &       +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0*
3274 <      &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
3275 <      &       sinphi/sinthet*costhet*costhet_grad(j)
3276 <      &      +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
3277 <      & )*wshield
3278 < C grad_shield_side is Cbeta sidechain gradient
3279 <       grad_shield_side(j,ishield_list(i),i)=
3280 <      &        (sh_frac_dist_grad(j)*-2.0d0
3281 <      &        *VofOverlap
3282 <      &       -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
3283 <      &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
3284 <      &       sinphi/sinthet*costhet*costhet_grad(j)
3285 <      &      +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
3286 <      &       )*wshield
3287 <
3288 <        grad_shield_loc(j,ishield_list(i),i)=
3289 <      &       scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
3290 <      &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(
3291 <      &       sinthet/sinphi*cosphi*cosphi_grad_loc(j)
3292 <      &        ))
3293 <      &        *wshield
3294 <       enddo
3295 <       VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
3296 <       enddo
3297 <       fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
3298 < C      write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
3299 <       enddo