2 < include 'COMMON.SPLITELE'
6 < C print *,"bylem w egb"
9 < cmc Sep-06: egb takes care of dynamic ss bonds too
11 < c if (dyn_ss) call dyn_set_nss
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
21 < else if (shield_mode.eq.2) then
22 < call set_shield_fac2
25 < write (iout,*) "Soft-spheer ELEC potential"
26 < c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
29 > c write (iout,*) "Soft-spheer ELEC potential"
30 > call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
33 < if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
34 < call ebend(ebe,ethetacnstr)
36 < C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
38 < if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
39 < call ebend_kcc(ebe,ethetacnstr)
46 < C print *,"TU DOCHODZE?"
48 < C print *,"tor",tor_mode
50 < if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
53 < C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
55 < if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
56 < call etor_kcc(etors,edihcnstr)
59 < if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then
61 > if (wtor_d.gt.0) then
63 < C print *,"PRZED MULIt"
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)
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)
79 < energia(22)=eliptran
80 < energia(23)=Eafmforce
81 < energia(24)=ethetacnstr
83 < if (dyn_ss) call dyn_set_nss
85 < eliptran=energia(22)
86 < Eafmforce=energia(23)
87 < ethetacnstr=energia(24)
89 < & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
92 > & +wbond*estr+Uconst+wsccor*esccor
94 < & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
98 > & +wbond*estr+Uconst+wsccor*esccor
100 > double precision gradbufc(3,maxres),gradbufx(3,maxres),
101 > & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
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)
111 < & +wliptran*gliptranc(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)
125 < & +wliptran*gliptranc(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)
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"
162 < & +wliptran*gliptranc(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)
181 < & welec*gelc_long(j,i)+
183 > & welec*gelc_long(j,i)
185 < & +wliptran*gliptranc(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)
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)
221 < include 'COMMON.CONTROL'
223 < if (shield_mode.gt.0) then
224 < wscp=weights(2)*fact
225 < wsc=weights(1)*fact
226 < wvdwpp=weights(16)*fact
229 < eliptran=energia(22)
230 < Eafmforce=energia(23)
231 < ethetacnstr=energia(24)
233 < & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
234 < & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
237 > & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
238 > & edihcnstr,ebr*nss,
241 < & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
243 < & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
244 < & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/
248 < & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
251 > & ebr*nss,Uconst,etot
253 < & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
255 < & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
256 < & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/
258 < C have you changed here?
262 > e1=fac*fac*aa(itypi,itypj)
263 > e2=fac*bb(itypi,itypj)
265 < cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj),
267 > cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
269 < C have you changed here?
273 > e1=fac*fac*aa(itypi,itypj)
274 > e2=fac*bb(itypi,itypj)
276 < C have you changed here?
281 > e1=fac*fac*aa(itypi,itypj)
282 > e2=fac*bb(itypi,itypj)
284 < sigm=dabs(aa/bb)**(1.0D0/6.0D0)
287 > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
288 > epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
290 < include 'COMMON.SPLITELE'
291 < include 'COMMON.SBRIDGE'
293 < integer xshift,yshift,zshift
296 < C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
298 > c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
300 < C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0
301 < C we have the original box)
306 < C Return atom into box, boxxsize is size of box in x dimension
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
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
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
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
339 < C if (positi.le.0) positi=positi+boxzsize
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
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
366 < C xi=xi+xshift*boxxsize
367 < C yi=yi+yshift*boxysize
368 < C zi=zi+zshift*boxzsize
371 < IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
373 < c write(iout,*) "PRZED ZWYKLE", evdwij
374 < call dyn_ssbond_ene(i,j,evdwij)
375 < c write(iout,*) "PO ZWYKLE", 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
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
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
398 < if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
399 < & 'evdw',i,j,evdwij,'tss'
400 < endif!dyn_ss_mask(k)
407 < C Return atom J into box the original box
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
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
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
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
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
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
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
490 < if (subchap.eq.1) then
508 < sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
509 < sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
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
515 < C here to start with
521 > e1=fac*fac*aa(itypi,itypj)
522 > e2=fac*bb(itypi,itypj)
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)
528 < evdw=evdw+evdwij*sss
532 < sigm=dabs(aa/bb)**(1.0D0/6.0D0)
535 > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
536 > epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
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
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
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
566 < C if (positi.le.0) positi=positi+boxzsize
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
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
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
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
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
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
653 < if (subchap.eq.1) then
670 > e1=fac*fac*aa(itypi,itypj)
671 > e2=fac*bb(itypi,itypj)
673 < sigm=dabs(aa/bb)**(1.0D0/6.0D0)
676 > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
677 > epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
679 < fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
681 < cc print *,'sss=',sss
683 < gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss
685 > gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
687 < gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
689 > gvdwx(k,i)=gvdwx(k,i)-gg(k)
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)
694 > & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
695 > gvdwx(k,j)=gvdwx(k,j)+gg(k)
697 < & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss
699 > & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
701 < gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
702 < gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
704 > gvdwc(l,i)=gvdwc(l,i)-gg(l)
705 > gvdwc(l,j)=gvdwc(l,j)+gg(l)
707 < C write(iout,*) 'In EELEC_soft_sphere'
709 > cd write(iout,*) 'In EELEC_soft_sphere'
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
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
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
749 < if (isubchap.eq.1) then
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
763 < sss=sscale(sqrt(rij))
764 < sssgrad=sscagrad(sqrt(rij))
766 < evdw1=evdw1+evdw1ij*sss
768 > evdw1=evdw1+evdw1ij
770 < ggg(1)=fac*xj*sssgrad
771 < ggg(2)=fac*yj*sssgrad
772 < ggg(3)=fac*zj*sssgrad
778 < c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
780 > c & +bnew1(3,1,iti)*dsin(alpha(i))*cos(beta(i))
782 < c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
784 > c & +bnew2(3,1,iti)*dsin(alpha(i))*dcos(beta(i))
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
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
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))
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))
800 < c write(iout,*) 'b1=',b1(1,i-2)
803 > do i=ivec_start+2,ivec_end+2
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)
824 < do i=ivec_start+2,ivec_end+2
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)
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)
835 < C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2)
836 < c write (iout,*) 'mu ',mu(:,i-2),i-2
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))
846 > cd write (iout,*) 'mu ',mu(:,i-2)
848 < cd iti = itortyp(itype(i))
850 > cd iti = itype2loc(itype(i))
852 < include 'COMMON.SPLITELE'
854 < C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
857 < C write(iout,*) "tu jest i",i
859 < C changes suggested by Ana to avoid out of bounds
860 < & .or.((i+4).gt.nres)
862 < C end of changes by Ana
863 < & .or. itype(i+2).eq.ntyp1
864 < & .or. itype(i+3).eq.ntyp1) cycle
866 < if(itype(i-1).eq.ntyp1)cycle
868 < if(i.LT.nres-3)then
869 < if (itype(i+4).eq.ntyp1) cycle
872 > & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
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
883 < C changes suggested by Ana to avoid out of bounds
884 < & .or.((i+5).gt.nres)
886 < C end of changes suggested by Ana
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
894 > & .or. itype(i+4).eq.ntyp1) cycle
896 < C Return atom into box, boxxsize is size of box in x dimension
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
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
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
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
929 < C Loop over all neighbouring boxes
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)
942 < C end of changes by Ana
943 < & .or. itype(i+2).eq.ntyp1
944 < & .or. itype(i-1).eq.ntyp1
948 > if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
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
960 < C Return tom into box, boxxsize is size of box in x dimension
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
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
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
990 < C write (iout,*) i,j
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)
996 < C end of changes by Ana
997 < & .or.itype(j+2).eq.ntyp1
998 < & .or.itype(j-1).eq.ntyp1
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
1010 < include 'COMMON.SPLITELE'
1011 < include 'COMMON.SHIELD'
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
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
1048 < if (isubchap.eq.1) then
1057 < C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
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
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
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
1082 < C endif !endPBC condintion
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
1092 < sss=sscale(sqrt(rij))
1093 < sssgrad=sscagrad(sqrt(rij))
1094 < c if (sss.gt.0.0d0) then
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
1118 < evdw1=evdw1+evdwij*sss
1120 > evdw1=evdw1+evdwij
1122 < write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
1123 < &fac_shield(i),fac_shield(j)
1125 > write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
1127 < facvdw=-6*rrmij*(ev1+evdwij)*sss
1129 > facvdw=-6*rrmij*(ev1+evdwij)
1133 < if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
1134 < & (shield_mode.gt.0)) then
1136 < do ilist=1,ishield_list(i)
1137 < iresshield=shield_list(ilist,i)
1139 < rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
1141 < gshieldx(k,iresshield)=gshieldx(k,iresshield)+
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)
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)
1162 < do ilist=1,ishield_list(j)
1163 < iresshield=shield_list(ilist,j)
1165 < rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
1167 < gshieldx(k,iresshield)=gshieldx(k,iresshield)+
1169 < & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
1170 < gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
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)
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)
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
1203 < C print *,"before", gelc_long(1,i), gelc_long(1,j)
1205 < C & +grad_shield(k,j)*eesij/fac_shield(j)
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)
1213 < C print *,"bafter", gelc_long(1,i), gelc_long(1,j)
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
1231 < facvdw=(ev1+evdwij)*sss
1237 < C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j)
1239 < C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j)
1241 < C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j)
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
1251 < ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
1252 < & fac_shield(i)**2*fac_shield(j)**2
1254 > ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
1256 < C print *,"before22", gelc_long(1,i), gelc_long(1,j)
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
1262 > & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
1263 > & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
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
1269 > & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
1270 > & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
1272 < C print *,"before33", gelc_long(1,i), gelc_long(1,j)
1277 < c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
1279 < c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
1281 > c write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
1283 < c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
1285 > c write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
1287 < cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
1289 > cd & i,itype2loc(itype(i)),j,itype2loc(itype(j)),a22,a23,a32,a33
1291 > c if ((i.eq.8).and.(j.eq.14)) then
1293 < if (shield_mode.eq.0) then
1297 < C fac_shield(i)=0.4
1298 < C fac_shield(j)=0.6
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
1307 < do ilist=1,ishield_list(i)
1308 < iresshield=shield_list(ilist,i)
1310 < rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
1313 < gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
1315 < & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
1316 < gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
1320 < do ilist=1,ishield_list(j)
1321 < iresshield=shield_list(ilist,j)
1323 < rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
1326 < gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
1328 < & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
1329 < gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
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)
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)
1352 < geel_loc_ij=(a22*gmuij1(1)
1354 > geel_loc_ij=a22*gmuij1(1)
1357 < & *fac_shield(i)*fac_shield(j)
1366 > geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
1368 < & *fac_shield(i)*fac_shield(j)
1370 < c Derivative over j residue
1371 < geel_loc_ji=a22*gmuji1(1)
1375 > geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
1377 < gloc(nphi+j,icg)=gloc(nphi+j,icg)+
1379 > gloc(nphi+j,icg)=gloc(nphi+j,icg)+
1381 < & *fac_shield(i)*fac_shield(j)
1388 > geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
1390 < & *fac_shield(i)*fac_shield(j)
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)
1396 > c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
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)
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)
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)
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)
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)
1417 > ggg(l)=agg(l,1)*muij(1)+
1418 > & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
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)
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)
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)
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)
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)
1448 < if (shield_mode.eq.0) then
1449 < fac_shield(i)=1.0d0
1450 < fac_shield(j)=1.0d0
1452 < ees0plist(num_conti,i)=j
1453 < C fac_shield(i)=0.4d0
1454 < C fac_shield(j)=0.6d0
1457 < & *fac_shield(i)*fac_shield(j)
1459 < & *fac_shield(i)*fac_shield(j)
1461 < & *fac_shield(i)*fac_shield(j)
1464 < & *fac_shield(i)*fac_shield(j)
1467 < & *fac_shield(i)*fac_shield(j)
1470 < & *fac_shield(i)*fac_shield(j)
1473 < & *fac_shield(i)*fac_shield(j)
1476 < & *fac_shield(i)*fac_shield(j)
1479 < include 'COMMON.SHIELD'
1481 < if (shield_mode.eq.0) then
1485 < C fac_shield(i)=0.4
1486 < C fac_shield(j)=0.6
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)
1493 < & *fac_shield(i)*fac_shield(j)
1495 < & *fac_shield(i)*fac_shield(j)
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
1504 < do ilist=1,ishield_list(i)
1505 < iresshield=shield_list(ilist,i)
1507 < rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i)
1509 < gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
1511 < & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i)
1512 < gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
1516 < do ilist=1,ishield_list(j)
1517 < iresshield=shield_list(ilist,j)
1519 < rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j)
1521 < gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
1523 < & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j)
1524 < gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
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)
1542 < C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1544 > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
1545 > & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
1547 < & *fac_shield(i)*fac_shield(j)
1549 < & *fac_shield(i)*fac_shield(j)
1551 < & *fac_shield(i)*fac_shield(j)
1554 < & *fac_shield(i)*fac_shield(j)
1556 < & *fac_shield(i)*fac_shield(j)
1558 < & *fac_shield(i)*fac_shield(j)
1560 < include 'COMMON.SHIELD'
1562 < if (shield_mode.eq.0) then
1566 < C fac_shield(i)=0.6
1567 < C fac_shield(j)=0.4
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
1581 < do ilist=1,ishield_list(i)
1582 < iresshield=shield_list(ilist,i)
1584 < rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
1586 < gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
1588 < & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)
1589 < gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
1593 < do ilist=1,ishield_list(j)
1594 < iresshield=shield_list(ilist,j)
1596 < rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j)
1598 < gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
1600 < & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j)
1601 < gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
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)
1624 < cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
1625 < cd & ' eello_turn4_num',8*eello_turn4_num
1627 > eello_turn4=eello_turn4-(s1+s2+s3)
1629 < & *fac_shield(i)*fac_shield(j)
1631 < & *fac_shield(i)*fac_shield(j)
1634 < & *fac_shield(i)*fac_shield(j)
1637 < & *fac_shield(i)*fac_shield(j)
1639 < & *fac_shield(i)*fac_shield(j)
1641 < & *fac_shield(i)*fac_shield(j)
1643 < & *fac_shield(i)*fac_shield(j)
1645 < & *fac_shield(i)*fac_shield(j)
1647 < & *fac_shield(i)*fac_shield(j)
1649 < & *fac_shield(i)*fac_shield(j)
1651 < & *fac_shield(i)*fac_shield(j)
1657 < C Return atom into box, boxxsize is size of box in x dimension
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
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
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
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
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
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
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
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
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
1748 < if (subchap.eq.1) then
1772 < include 'COMMON.SPLITELE'
1774 < c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
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
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
1800 < c print *,xi,boxxsize,"pierwszy"
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
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
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
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
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
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
1851 < CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
1852 < dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
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
1874 < if (subchap.eq.1) then
1883 < c print *,xj,yj,zj,'polozenie j'
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)))
1896 < evdw2_14=evdw2_14+(e1+e2)*sss
1898 > evdw2_14=evdw2_14+e1+e2
1900 < evdw2=evdw2+evdwij*sss
1902 > evdw2=evdw2+evdwij
1904 < fac=-(evdwij+e1)*rrij*sss
1905 < fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
1907 > fac=-(evdwij+e1)*rrij
1909 < c endif !endif for sscale cutoff
1918 < include 'COMMON.CONTROL'
1923 < C write (iout,*) ,"link_end",link_end,constr_dist
1925 < c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
1926 < c & dhpb(i),dhpb1(i),forcon(i)
1928 > cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
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.
1938 > if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
1942 < cd & ' waga=',waga,' fac=',fac
1943 < else if (ii.gt.nres .and. jj.gt.nres) then
1944 < c Restraints from contact prediction
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
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))
1962 < C Get the force constant corresponding to this distance.
1964 < C Calculate the contribution to energy.
1965 < ehpb=ehpb+waga*rdis*rdis
1966 < c write (iout,*) "beta reg",dd,waga*rdis*rdis
1968 < C Evaluate gradient.
1974 < ggg(j)=fac*(c(j,jj)-c(j,ii))
1977 < ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
1978 < ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
1981 < ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
1982 < ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
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
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))
2009 < ehpb=ehpb+waga*rdis*rdis
2010 < c write (iout,*) "alpha reg",dd,waga*rdis*rdis
2012 > ehpb=ehpb+waga*rdis*rdis
2018 < ggg(j)=fac*(c(j,jj)-c(j,ii))
2022 > cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
2023 > cd & ' waga=',waga,' fac=',fac
2025 > ggg(j)=fac*(c(j,jj)-c(j,ii))
2028 < if (iii.lt.ii) then
2030 > if (iii.lt.ii) then
2037 < ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
2038 < ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
2042 > ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
2043 > ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
2046 < if (constr_dist.ne.11) ehpb=0.5D0*ehpb
2050 < if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
2051 < c estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
2053 < c gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
2054 < c & *dc(j,i-1)/vbld(i)
2056 < c if (energy_dec) write(iout,*)
2057 < c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
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
2064 < C NO vbldp0 is the equlibrium lenght of spring for peptide group
2066 > if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
2067 > estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
2069 > gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
2070 > & *dc(j,i-1)/vbld(i)
2072 > if (energy_dec) write(iout,*)
2073 > & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
2077 < if (energy_dec) write (iout,'(a7,i5,4f7.3)')
2079 > if (energy_dec) write (iout,*)
2085 < if (energy_dec) write (iout,*)
2087 > if (energy_dec) write (iout,*)
2089 < subroutine ebend(etheta,ethetacnstr)
2091 > subroutine ebend(etheta)
2093 < include 'COMMON.TORCNSTR'
2095 < if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
2096 < & .or.itype(i).eq.ntyp1) cycle
2098 > if (itype(i-1).eq.ntyp1) cycle
2100 < if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
2102 > if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
2104 < if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
2106 > if (i.lt.nres .and. itype(i).ne.ntyp1) then
2112 < c write(iout,*) 'chuj tu', y(k),z(k)
2114 < if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
2115 < & 'ebend',i,ethetai,theta(i),itype(i)
2117 > if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
2118 > & 'ebend',i,ethetai
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
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)
2149 < C the distributioni.
2150 < ccc write (iout,*) thetai,thet_pred_mean
2152 > C the distribution.
2154 < C write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0
2156 < C write (iout,*) 'termexp',termexp,termm,termpre,i
2158 < subroutine ebend(etheta,ethetacnstr)
2160 > subroutine ebend(etheta)
2162 < include 'COMMON.TORCNSTR'
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)
2169 > if (itype(i-1).eq.ntyp1) cycle
2172 < if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
2174 > if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
2178 < ityp1=ithetyp((itype(i-2)))
2180 < if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
2182 > if (i.lt.nres .and. itype(i).ne.ntyp1) then
2184 < ityp3=ithetyp((itype(i)))
2190 < C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
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)
2197 < C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
2201 < gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
2203 > gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
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
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)
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))
2239 > & .or. itype(i).eq.ntyp1) cycle
2240 > itori=itortyp(itype(i-2))
2241 > itori1=itortyp(itype(i-1))
2243 < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2244 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2246 > edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2247 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
2249 < edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4
2250 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2252 > edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2253 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
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
2265 > if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
2266 > & .or. itype(i).eq.ntyp1) cycle
2268 < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2269 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2271 > edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2272 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
2274 < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
2275 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
2277 > edihcnstr=edihcnstr+0.25d0*ftors*difi**4
2278 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
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)
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)
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
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
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
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)
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'
2336 < c double precision thybt1(maxtermkcc),thybt2(maxtermkcc)
2337 < C Set lprn=.true. for debugging
2340 < C print *,"wchodze kcc"
2341 < if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode
2342 < if (tor_mode.ne.2) then
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))
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
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
2381 < do j=1,nterm_kcc(itori,itori1)
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),
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),
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),
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),
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)
2424 < & j*sint1t2n*(v2ij*cosphi*(1.0d0+sumth1tyb2+sumth2tyb2)
2425 < & -v1ij*sinphi*(1.0d0+sumth1tyb1+sumth2tyb1))
2426 < C now gradient over theta_1
2428 < & j*sint1t2n1*costhet1*sinthet2*(actval1+actval2)+
2429 < & sint1t2n*(v1ij*cosphi*gradth1tyb1+v2ij*sinphi*gradth1tyb2)
2431 < & j*sint1t2n1*sinthet1*costhet2*(actval1+actval2)+
2432 < & sint1t2n*(v1ij*cosphi*gradth2tyb1+v2ij*sinphi*gradth2tyb2)
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)
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),
2447 < C now overal sumation
2448 < C print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb
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
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
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
2465 < c do i=1,ndih_constr
2466 < do i=idihconstr_start,idihconstr_end
2467 < itori=idih_constr(i)
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
2486 < C The rigorous attempt to derive energy function
2487 < subroutine ebend_kcc(etheta,ethetacnstr)
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'
2504 < double precision thybt1(maxtermkcc)
2505 < C Set lprn=.true. for debugging
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)
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,
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)
2533 < if (tor_mode.ne.2) then
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
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)
2565 < include 'COMMON.SHIELD'
2567 < include 'COMMON.SHIELD'
2569 < CC & *fac_shield(i)**2*fac_shield(j)**2
2571 < include 'COMMON.SHIELD'
2572 < include 'COMMON.CONTROL'
2574 < C print *,"wchodze",fac_shield(i),shield_mode
2577 < C & fac_shield(i)**2*fac_shield(j)**2
2579 < C print *,ekont,ees,i,k
2581 < C now gradient over shielding
2583 < if (shield_mode.gt.0) then
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)
2593 < rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i)
2595 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2597 < & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i)
2598 < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2602 < do ilist=1,ishield_list(j)
2603 < iresshield=shield_list(ilist,j)
2605 < rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j)
2607 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2609 < & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j)
2610 < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2615 < do ilist=1,ishield_list(k)
2616 < iresshield=shield_list(ilist,k)
2618 < rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k)
2620 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2622 < & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k)
2623 < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2627 < do ilist=1,ishield_list(l)
2628 < iresshield=shield_list(ilist,l)
2630 < rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l)
2632 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
2634 < & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l)
2635 < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
2639 < C print *,gshieldx(m,iresshield)
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)
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)
2663 < iti1 = itortyp(itype(i+1))
2665 > iti1 = itype2loc(itype(i+1))
2667 < call transpose2(EE(1,1,k),auxmat(1,1))
2669 > call transpose2(EE(1,1,itk),auxmat(1,1))
2671 < call transpose2(EE(1,1,l),auxmat(1,1))
2673 > call transpose2(EE(1,1,itl),auxmat(1,1))
2675 < call transpose2(EE(1,1,j),auxmat(1,1))
2677 > call transpose2(EE(1,1,itj),auxmat(1,1))
2679 < gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2)
2681 > gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2)
2683 < gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2)
2685 > gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2)
2690 10084,10085c8253,8254
2705 < iti=itortyp(itype(i))
2707 > iti=itype2loc(itype(i))
2709 < call transpose2(EE(1,1,k),auxmat(1,1))
2711 > call transpose2(EE(1,1,itk),auxmat(1,1))
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:
2747 < C--bordliptop-- buffore starts
2748 < C--bufliptop--- here true lipid starts
2750 < C--buflipbot--- lipid ends buffore starts
2751 < C--bordlipbot--buffore ends
2753 < do i=ilip_start,ilip_end
2755 < if (itype(i).eq.ntyp1) cycle
2757 < positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
2758 < if (positi.le.0) positi=positi+boxzsize
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
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
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
2790 < eliptran=eliptran+pepliptran
2791 < C print *,"I am in true lipid"
2794 < C eliptran=elpitran+0.0 ! I am in water
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
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
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
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
2837 < eliptran=eliptran+liptranene(itype(i))
2838 < C print *,"I am in true lipid"
2840 < endif ! if in lipid or buffor
2842 < C eliptran=elpitran+0.0 ! I am in water
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'
2867 < diffafm(i)=c(i,afmend)-c(i,afmbeg)
2868 < dist=dist+diffafm(i)**2
2871 < Eafmforce=-forceAFMconst*(dist-distafminit)
2873 < gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist
2874 < gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist
2876 < C print *,'AFM',Eafmforce
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'
2897 < C Only for check grad COMMENT if not used for checkgrad
2899 < C--------------------------------------------------------
2900 < C print *,"wchodze"
2904 < diffafm(i)=c(i,afmend)-c(i,afmbeg)
2905 < dist=dist+diffafm(i)**2
2908 < Eafmforce=0.5d0*forceAFMconst
2909 < & *(distafminit+totTafm*velAFMconst-dist)**2
2910 < C Eafmforce=-forceAFMconst*(dist-distafminit)
2912 < gradafm(i,afmend-1)=-forceAFMconst*
2913 < &(distafminit+totTafm*velAFMconst-dist)
2915 < gradafm(i,afmbeg-1)=forceAFMconst*
2916 < &(distafminit+totTafm*velAFMconst-dist)
2919 < C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
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)
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
2942 < if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
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
2948 < if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
2950 < dist_side_calf=0.0
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
2964 < dist_pep_side=dsqrt(dist_pep_side)
2965 < dist_pept_group=dsqrt(dist_pept_group)
2966 < dist_side_calf=dsqrt(dist_side_calf)
2968 < pep_side_norm(j)=pep_side(j)/dist_pep_side
2969 < side_calf_norm(j)=dist_side_calf
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"
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
2986 < sh_frac_dist_grad(j)=0.0d0
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 !
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)
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
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
3012 < costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
3013 < C costhet_fac=0.0d0
3015 < costhet_grad(j)=costhet_fac*pep_side(j)
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
3023 < pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
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
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))
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))*
3046 < &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
3049 < VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
3050 < & /VSolvSphere_div
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)
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
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
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
3080 < VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
3082 < fac_shield(i)=VolumeTotal*div77_81+div4_81
3083 < C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
3087 < C--------------------------------------------------------------------------
3088 < double precision function tschebyshev(m,n,x,y)
3090 < include "DIMENSIONS"
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
3099 < yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
3103 < aux=aux+x(i)*yy(i)
3108 < C--------------------------------------------------------------------------
3109 < double precision function gradtschebyshev(m,n,x,y)
3111 < include "DIMENSIONS"
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
3120 < yy(i)=2*y*yy(i-1)-yy(i-2)
3124 < aux=aux+x(i+1)*yy(i)*(i+1)
3125 < C print *, x(i+1),yy(i),i
3127 < gradtschebyshev=aux
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)
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
3150 < if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
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
3156 < if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
3158 < dist_side_calf=0.0
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
3172 < dist_pep_side=dsqrt(dist_pep_side)
3173 < dist_pept_group=dsqrt(dist_pept_group)
3174 < dist_side_calf=dsqrt(dist_side_calf)
3176 < pep_side_norm(j)=pep_side(j)/dist_pep_side
3177 < side_calf_norm(j)=dist_side_calf
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"
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
3194 < sh_frac_dist_grad(j)=0.0d0
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 !
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)
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
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
3224 < costhet_grad(j)=costhet_fac*pep_side(j)
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
3232 < pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
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
3241 < C now costhet_grad
3242 < cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
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)
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))*
3260 < &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
3261 < C cosphi_grad_loc(j)=0.0d0
3263 < C print *,sinphi,sinthet
3264 < VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
3265 < & /VSolvSphere_div
3267 < C now the gradient...
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)))
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
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)))
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)
3295 < VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
3297 < fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
3298 < C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)