27d26 < include 'COMMON.SPLITELE' 102d100 < C print *,ipot 116d113 < C print *,"bylem w egb" 127,131d123 < cmc < cmc Sep-06: egb takes care of dynamic ss bonds too < cmc < c if (dyn_ss) call dyn_set_nss < 140,149d131 < C Introduction of shielding effect first for each peptide group < C the shielding factor is set this factor is describing how each < C peptide group is shielded by side-chains < C the matrix - shield_fac(i) the i index describe the ith between i and i+1 < C write (iout,*) "shield_mode",shield_mode < if (shield_mode.eq.1) then < call set_shield_fac < else if (shield_mode.eq.2) then < call set_shield_fac2 < endif 172,174c154,156 < write (iout,*) "Soft-spheer ELEC potential" < c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, < c & eello_turn4) --- > c write (iout,*) "Soft-spheer ELEC potential" > call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, > & eello_turn4) 206,213c188 < if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then < call ebend(ebe,ethetacnstr) < endif < C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the < C energy function < if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then < call ebend_kcc(ebe,ethetacnstr) < endif --- > call ebend(ebe) 216d190 < ethetacnstr=0 222d195 < C print *,"TU DOCHODZE?" 229d201 < C print *,"tor",tor_mode 231d202 < if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then 233,238d203 < endif < C etor kcc is Kubo cumulant clustered rigorous attemp to derive the < C energy function < if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then < call etor_kcc(etors,edihcnstr) < endif 247c212 < if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then --- > if (wtor_d.gt.0) then 261d225 < C print *,"PRZED MULIt" 294,306d257 < C 01/27/2015 added by adasko < C the energy component below is energy transfer into lipid environment < C based on partition function < C print *,"przed lipidami" < if (wliptran.gt.0) then < call Eliptransfer(eliptran) < endif < C print *,"za lipidami" < if (AFMlog.gt.0) then < call AFMforce(Eafmforce) < else if (selfguide.gt.0) then < call AFMvel(Eafmforce) < endif 348,350d298 < energia(22)=eliptran < energia(23)=Eafmforce < energia(24)=ethetacnstr 355d302 < if (dyn_ss) call dyn_set_nss 442,444d388 < eliptran=energia(22) < Eafmforce=energia(23) < ethetacnstr=energia(24) 451,452c395 < & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce < & +ethetacnstr --- > & +wbond*estr+Uconst+wsccor*esccor 459,461c402 < & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran < & +Eafmforce < & +ethetacnstr --- > & +wbond*estr+Uconst+wsccor*esccor 496a438,439 > double precision gradbufc(3,maxres),gradbufx(3,maxres), > & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) 498,500d440 < double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres), < & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres) < & ,gloc_scbuf(3,-1:maxres) 553c493 < do i=0,nct --- > do i=1,nct 564,572d503 < & +wliptran*gliptranc(j,i) < & +gradafm(j,i) < & +welec*gshieldc(j,i) < & +wcorr*gshieldc_ec(j,i) < & +wturn3*gshieldc_t3(j,i) < & +wturn4*gshieldc_t4(j,i) < & +wel_loc*gshieldc_ll(j,i) < < 576c507 < do i=0,nct --- > do i=1,nct 588,595d518 < & +wliptran*gliptranc(j,i) < & +gradafm(j,i) < & +welec*gshieldc(j,i) < & +wcorr*gshieldc_ec(j,i) < & +wturn4*gshieldc_t4(j,i) < & +wel_loc*gshieldc_ll(j,i) < < 609c532 < do i=0,nres --- > do i=1,nres 652c575 < do i=nres-2,-1,-1 --- > do i=nres-2,nnt,-1 673c596 < do i=-1,nres --- > do i=1,nres 682c605 < do i=nres-2,-1,-1 --- > do i=nres-2,nnt,-1 710c633 < do i=-1,nct --- > do i=1,nct 713,719d635 < C print *,gradbufc(1,13) < C print *,welec*gelc(1,13) < C print *,wel_loc*gel_loc(1,13) < C print *,0.5d0*(wscp*gvdwc_scpp(1,13)) < C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13) < C print *,wel_loc*gel_loc_long(1,13) < C print *,gradafm(1,13),"AFM" 738,755d653 < & +wliptran*gliptranc(j,i) < & +gradafm(j,i) < & +welec*gshieldc(j,i) < & +welec*gshieldc_loc(j,i) < & +wcorr*gshieldc_ec(j,i) < & +wcorr*gshieldc_loc_ec(j,i) < & +wturn3*gshieldc_t3(j,i) < & +wturn3*gshieldc_loc_t3(j,i) < & +wturn4*gshieldc_t4(j,i) < & +wturn4*gshieldc_loc_t4(j,i) < & +wel_loc*gshieldc_ll(j,i) < & +wel_loc*gshieldc_loc_ll(j,i) < < < < < < 760c658 < & welec*gelc_long(j,i)+ --- > & welec*gelc_long(j,i) 775,791d672 < & +wliptran*gliptranc(j,i) < & +gradafm(j,i) < & +welec*gshieldc(j,i) < & +welec*gshieldc_loc(j,i) < & +wcorr*gshieldc_ec(j,i) < & +wcorr*gshieldc_loc_ec(j,i) < & +wturn3*gshieldc_t3(j,i) < & +wturn3*gshieldc_loc_t3(j,i) < & +wturn4*gshieldc_t4(j,i) < & +wturn4*gshieldc_loc_t4(j,i) < & +wel_loc*gshieldc_ll(j,i) < & +wel_loc*gshieldc_loc_ll(j,i) < < < < < 798,806d678 < & +wliptran*gliptranx(j,i) < & +welec*gshieldx(j,i) < & +wcorr*gshieldx_ec(j,i) < & +wturn3*gshieldx_t3(j,i) < & +wturn4*gshieldx_t4(j,i) < & +wel_loc*gshieldx_ll(j,i) < < < 841d712 < c#define DEBUG 850d720 < c#undef DEBUG 870d739 < c#define DEBUG 879d747 < c#undef DEBUG 1009d876 < include 'COMMON.CONTROL' 1045,1049d911 < if (shield_mode.gt.0) then < wscp=weights(2)*fact < wsc=weights(1)*fact < wvdwpp=weights(16)*fact < endif 1101,1103d962 < eliptran=energia(22) < Eafmforce=energia(23) < ethetacnstr=energia(24) 1110,1112c969,971 < & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, < & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, < & etot --- > & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, > & edihcnstr,ebr*nss, > & Uconst,etot 1134d992 < & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ 1137,1138d994 < & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ < & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ 1140d995 < 1148,1149c1003 < & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, < & etot --- > & ebr*nss,Uconst,etot 1170d1023 < & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ 1173,1174d1025 < & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ < & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ 1229,1231c1080,1081 < C have you changed here? < e1=fac*fac*aa < e2=fac*bb --- > e1=fac*fac*aa(itypi,itypj) > e2=fac*bb(itypi,itypj) 1236c1086 < cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj), --- > cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), 1380,1382c1230,1231 < C have you changed here? < e1=fac*fac*aa < e2=fac*bb --- > e1=fac*fac*aa(itypi,itypj) > e2=fac*bb(itypi,itypj) 1508d1356 < C have you changed here? 1510,1511c1358,1359 < e1=fac*fac*aa < e2=fac*bb --- > e1=fac*fac*aa(itypi,itypj) > e2=fac*bb(itypi,itypj) 1518,1519c1366,1367 < sigm=dabs(aa/bb)**(1.0D0/6.0D0) < epsi=bb**2/aa --- > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) > epsi=bb(itypi,itypj)**2/aa(itypi,itypj) 1563,1564d1410 < include 'COMMON.SPLITELE' < include 'COMMON.SBRIDGE' 1566,1567d1411 < integer xshift,yshift,zshift < 1570c1414 < C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon --- > c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon 1575,1579d1418 < C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0 < C we have the original box) < C do xshift=-1,1 < C do yshift=-1,1 < C do zshift=-1,1 1587,1650d1425 < C Return atom into box, boxxsize is size of box in x dimension < c 134 continue < c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize < c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize < C Condition for being inside the proper box < c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. < c & (xi.lt.((xshift-0.5d0)*boxxsize))) then < c go to 134 < c endif < c 135 continue < c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize < c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize < C Condition for being inside the proper box < c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. < c & (yi.lt.((yshift-0.5d0)*boxysize))) then < c go to 135 < c endif < c 136 continue < c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize < c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize < C Condition for being inside the proper box < c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. < c & (zi.lt.((zshift-0.5d0)*boxzsize))) then < c go to 136 < c endif < xi=mod(xi,boxxsize) < if (xi.lt.0) xi=xi+boxxsize < yi=mod(yi,boxysize) < if (yi.lt.0) yi=yi+boxysize < zi=mod(zi,boxzsize) < if (zi.lt.0) zi=zi+boxzsize < C define scaling factor for lipids < < C if (positi.le.0) positi=positi+boxzsize < C print *,i < C first for peptide groups < c for each residue check if it is in lipid or lipid water border area < if ((zi.gt.bordlipbot) < &.and.(zi.lt.bordliptop)) then < C the energy transfer exist < if (zi.lt.buflipbot) then < C what fraction I am in < fracinbuf=1.0d0- < & ((zi-bordlipbot)/lipbufthick) < C lipbufthick is thickenes of lipid buffore < sslipi=sscalelip(fracinbuf) < ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick < elseif (zi.gt.bufliptop) then < fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) < sslipi=sscalelip(fracinbuf) < ssgradlipi=sscagradlip(fracinbuf)/lipbufthick < else < sslipi=1.0d0 < ssgradlipi=0.0 < endif < else < sslipi=0.0d0 < ssgradlipi=0.0 < endif < < C xi=xi+xshift*boxxsize < C yi=yi+yshift*boxysize < C zi=zi+zshift*boxzsize < 1663,1694d1437 < IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN < < c write(iout,*) "PRZED ZWYKLE", evdwij < call dyn_ssbond_ene(i,j,evdwij) < c write(iout,*) "PO ZWYKLE", evdwij < < evdw=evdw+evdwij < if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') < & 'evdw',i,j,evdwij,' ss' < C triple bond artifac removal < do k=j+1,iend(i,iint) < C search over all next residues < if (dyn_ss_mask(k)) then < C check if they are cysteins < C write(iout,*) 'k=',k < < c write(iout,*) "PRZED TRI", evdwij < evdwij_przed_tri=evdwij < call triple_ssbond_ene(i,j,k,evdwij) < c if(evdwij_przed_tri.ne.evdwij) then < c write (iout,*) "TRI:", evdwij, evdwij_przed_tri < c endif < < c write(iout,*) "PO TRI", evdwij < C call the energy function that removes the artifical triple disulfide < C bond the soubroutine is located in ssMD.F < evdw=evdw+evdwij < if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') < & 'evdw',i,j,evdwij,'tss' < endif!dyn_ss_mask(k) < enddo! k < ELSE 1723,1817c1466,1468 < xj=c(1,nres+j) < yj=c(2,nres+j) < zj=c(3,nres+j) < C Return atom J into box the original box < c 137 continue < c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize < c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize < C Condition for being inside the proper box < c if ((xj.gt.((0.5d0)*boxxsize)).or. < c & (xj.lt.((-0.5d0)*boxxsize))) then < c go to 137 < c endif < c 138 continue < c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize < c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize < C Condition for being inside the proper box < c if ((yj.gt.((0.5d0)*boxysize)).or. < c & (yj.lt.((-0.5d0)*boxysize))) then < c go to 138 < c endif < c 139 continue < c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize < c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize < C Condition for being inside the proper box < c if ((zj.gt.((0.5d0)*boxzsize)).or. < c & (zj.lt.((-0.5d0)*boxzsize))) then < c go to 139 < c endif < xj=mod(xj,boxxsize) < if (xj.lt.0) xj=xj+boxxsize < yj=mod(yj,boxysize) < if (yj.lt.0) yj=yj+boxysize < zj=mod(zj,boxzsize) < if (zj.lt.0) zj=zj+boxzsize < if ((zj.gt.bordlipbot) < &.and.(zj.lt.bordliptop)) then < C the energy transfer exist < if (zj.lt.buflipbot) then < C what fraction I am in < fracinbuf=1.0d0- < & ((zj-bordlipbot)/lipbufthick) < C lipbufthick is thickenes of lipid buffore < sslipj=sscalelip(fracinbuf) < ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick < elseif (zj.gt.bufliptop) then < fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) < sslipj=sscalelip(fracinbuf) < ssgradlipj=sscagradlip(fracinbuf)/lipbufthick < else < sslipj=1.0d0 < ssgradlipj=0.0 < endif < else < sslipj=0.0d0 < ssgradlipj=0.0 < endif < aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 < & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 < bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 < & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 < C if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)') < C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) < C if (ssgradlipj.gt.0.0d0) print *,"??WTF??" < C print *,sslipi,sslipj,bordlipbot,zi,zj < dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < xj_safe=xj < yj_safe=yj < zj_safe=zj < subchap=0 < do xshift=-1,1 < do yshift=-1,1 < do zshift=-1,1 < xj=xj_safe+xshift*boxxsize < yj=yj_safe+yshift*boxysize < zj=zj_safe+zshift*boxzsize < dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < if(dist_temp.lt.dist_init) then < dist_init=dist_temp < xj_temp=xj < yj_temp=yj < zj_temp=zj < subchap=1 < endif < enddo < enddo < enddo < if (subchap.eq.1) then < xj=xj_temp-xi < yj=yj_temp-yi < zj=zj_temp-zi < else < xj=xj_safe-xi < yj=yj_safe-yi < zj=zj_safe-zi < endif --- > xj=c(1,nres+j)-xi > yj=c(2,nres+j)-yi > zj=c(3,nres+j)-zi 1821,1823d1471 < C xj=xj-xi < C yj=yj-yi < C zj=zj-zi 1829,1834d1476 < sss=sscale((1.0d0/rij)/sigma(itypi,itypj)) < sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) < < c write (iout,'(a7,4f8.3)') < c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb < if (sss.gt.0.0d0) then 1855,1859c1497,1498 < C here to start with < C if (c(i,3).gt. < faclip=fac < e1=fac*fac*aa < e2=fac*bb --- > e1=fac*fac*aa(itypi,itypj) > e2=fac*bb(itypi,itypj) 1863,1865d1501 < C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij, < C &((sslipi+sslipj)/2.0d0+ < C &(2.0d0-sslipi-sslipj)/2.0d0) 1869c1505 < evdw=evdw+evdwij*sss --- > evdw=evdw+evdwij 1871,1872c1507,1508 < sigm=dabs(aa/bb)**(1.0D0/6.0D0) < epsi=bb**2/aa --- > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) > epsi=bb(itypi,itypj)**2/aa(itypi,itypj) 1889,1891d1524 < c print '(2i4,6f8.4)',i,j,sss,sssgrad* < c & evdwij,fac,sigma(itypi,itypj),expon < fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij 1894,1901d1526 < gg_lipi(3)=eps1*(eps2rt*eps2rt) < &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* < & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) < &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) < gg_lipj(3)=ssgradlipj*gg_lipi(3) < gg_lipi(3)=gg_lipi(3)*ssgradlipi < C gg_lipi(3)=0.0d0 < C gg_lipj(3)=0.0d0 1907,1908d1531 < endif < ENDIF ! dyn_ss 1912,1914d1534 < C enddo ! zshift < C enddo ! yshift < C enddo ! xshift 1951,1985d1570 < xi=mod(xi,boxxsize) < if (xi.lt.0) xi=xi+boxxsize < yi=mod(yi,boxysize) < if (yi.lt.0) yi=yi+boxysize < zi=mod(zi,boxzsize) < if (zi.lt.0) zi=zi+boxzsize < C define scaling factor for lipids < < C if (positi.le.0) positi=positi+boxzsize < C print *,i < C first for peptide groups < c for each residue check if it is in lipid or lipid water border area < if ((zi.gt.bordlipbot) < &.and.(zi.lt.bordliptop)) then < C the energy transfer exist < if (zi.lt.buflipbot) then < C what fraction I am in < fracinbuf=1.0d0- < & ((zi-bordlipbot)/lipbufthick) < C lipbufthick is thickenes of lipid buffore < sslipi=sscalelip(fracinbuf) < ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick < elseif (zi.gt.bufliptop) then < fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) < sslipi=sscalelip(fracinbuf) < ssgradlipi=sscagradlip(fracinbuf)/lipbufthick < else < sslipi=1.0d0 < ssgradlipi=0.0 < endif < else < sslipi=0.0d0 < ssgradlipi=0.0 < endif < 2022,2089c1607,1609 < C xj=c(1,nres+j)-xi < C yj=c(2,nres+j)-yi < C zj=c(3,nres+j)-zi < xj=mod(xj,boxxsize) < if (xj.lt.0) xj=xj+boxxsize < yj=mod(yj,boxysize) < if (yj.lt.0) yj=yj+boxysize < zj=mod(zj,boxzsize) < if (zj.lt.0) zj=zj+boxzsize < if ((zj.gt.bordlipbot) < &.and.(zj.lt.bordliptop)) then < C the energy transfer exist < if (zj.lt.buflipbot) then < C what fraction I am in < fracinbuf=1.0d0- < & ((zj-bordlipbot)/lipbufthick) < C lipbufthick is thickenes of lipid buffore < sslipj=sscalelip(fracinbuf) < ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick < elseif (zj.gt.bufliptop) then < fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) < sslipj=sscalelip(fracinbuf) < ssgradlipj=sscagradlip(fracinbuf)/lipbufthick < else < sslipj=1.0d0 < ssgradlipj=0.0 < endif < else < sslipj=0.0d0 < ssgradlipj=0.0 < endif < aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 < & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 < bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 < & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 < C if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5') < C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) < dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < xj_safe=xj < yj_safe=yj < zj_safe=zj < subchap=0 < do xshift=-1,1 < do yshift=-1,1 < do zshift=-1,1 < xj=xj_safe+xshift*boxxsize < yj=yj_safe+yshift*boxysize < zj=zj_safe+zshift*boxzsize < dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < if(dist_temp.lt.dist_init) then < dist_init=dist_temp < xj_temp=xj < yj_temp=yj < zj_temp=zj < subchap=1 < endif < enddo < enddo < enddo < if (subchap.eq.1) then < xj=xj_temp-xi < yj=yj_temp-yi < zj=zj_temp-zi < else < xj=xj_safe-xi < yj=yj_safe-yi < zj=zj_safe-zi < endif --- > xj=c(1,nres+j)-xi > yj=c(2,nres+j)-yi > zj=c(3,nres+j)-zi 2110,2111c1630,1631 < e1=fac*fac*aa < e2=fac*bb --- > e1=fac*fac*aa(itypi,itypj) > e2=fac*bb(itypi,itypj) 2120,2121c1640,1641 < sigm=dabs(aa/bb)**(1.0D0/6.0D0) < epsi=bb**2/aa --- > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) > epsi=bb(itypi,itypj)**2/aa(itypi,itypj) 2135d1654 < fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij 2223d1741 < cc print *,'sss=',sss 2242c1760 < gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss --- > gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k) 2246c1764 < gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) --- > gvdwx(k,i)=gvdwx(k,i)-gg(k) 2248,2249c1766,1767 < & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss < gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) --- > & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv > gvdwx(k,j)=gvdwx(k,j)+gg(k) 2251c1769 < & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss --- > & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv 2266,2267c1784,1785 < gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) < gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) --- > gvdwc(l,i)=gvdwc(l,i)-gg(l) > gvdwc(l,j)=gvdwc(l,j)+gg(l) 2369c1887 < C write(iout,*) 'In EELEC_soft_sphere' --- > cd write(iout,*) 'In EELEC_soft_sphere' 2384,2389d1901 < xmedi=mod(xmedi,boxxsize) < if (xmedi.lt.0) xmedi=xmedi+boxxsize < ymedi=mod(ymedi,boxysize) < if (ymedi.lt.0) ymedi=ymedi+boxysize < zmedi=mod(zmedi,boxzsize) < if (zmedi.lt.0) zmedi=zmedi+boxzsize 2403,2442c1915,1917 < xj=c(1,j)+0.5D0*dxj < yj=c(2,j)+0.5D0*dyj < zj=c(3,j)+0.5D0*dzj < xj=mod(xj,boxxsize) < if (xj.lt.0) xj=xj+boxxsize < yj=mod(yj,boxysize) < if (yj.lt.0) yj=yj+boxysize < zj=mod(zj,boxzsize) < if (zj.lt.0) zj=zj+boxzsize < dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 < xj_safe=xj < yj_safe=yj < zj_safe=zj < isubchap=0 < do xshift=-1,1 < do yshift=-1,1 < do zshift=-1,1 < xj=xj_safe+xshift*boxxsize < yj=yj_safe+yshift*boxysize < zj=zj_safe+zshift*boxzsize < dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < if(dist_temp.lt.dist_init) then < dist_init=dist_temp < xj_temp=xj < yj_temp=yj < zj_temp=zj < isubchap=1 < endif < enddo < enddo < enddo < if (isubchap.eq.1) then < xj=xj_temp-xmedi < yj=yj_temp-ymedi < zj=zj_temp-zmedi < else < xj=xj_safe-xmedi < yj=yj_safe-ymedi < zj=zj_safe-zmedi < endif --- > xj=c(1,j)+0.5D0*dxj-xmedi > yj=c(2,j)+0.5D0*dyj-ymedi > zj=c(3,j)+0.5D0*dzj-zmedi 2444,2445d1918 < sss=sscale(sqrt(rij)) < sssgrad=sscagrad(sqrt(rij)) 2453c1926 < evdw1=evdw1+evdw1ij*sss --- > evdw1=evdw1+evdw1ij 2457,2459c1930,1932 < ggg(1)=fac*xj*sssgrad < ggg(2)=fac*yj*sssgrad < ggg(3)=fac*zj*sssgrad --- > ggg(1)=fac*xj > ggg(2)=fac*yj > ggg(3)=fac*zj 2803c2276 < c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) --- > c & +bnew1(3,1,iti)*dsin(alpha(i))*cos(beta(i)) 2808c2281 < c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) --- > c & +bnew2(3,1,iti)*dsin(alpha(i))*dcos(beta(i)) 2815,2817c2288,2290 < c &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0, < c &bnew1(2,1,iti)*cos(theta(i)), < c &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0 --- > c &bnew1(1,1,iti)*dcos(theta(i)/2.0d0)/2.0d0, > c &bnew1(2,1,iti)*dcos(theta(i)), > c &bnew1(3,1,iti)*dsin(theta(i)/2.0d0)/2.0d0 2834,2835c2307,2308 < c b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i)) < c b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i)) --- > c b1(2,iti)=bnew1(1,2,iti)*dsin(alpha(i))*dsin(beta(i)) > c b2(2,iti)=bnew2(1,2,iti)*dsin(alpha(i))*dsin(beta(i)) 2841d2313 < c write(iout,*) 'b1=',b1(1,i-2) 2843a2316,2317 > #ifdef PARMAT > do i=ivec_start+2,ivec_end+2 2844a2319,2321 > do i=3,nres+1 > #endif > #endif 2856,2874d2332 < b1(1,i-2)=b(3,iti) < b1(2,i-2)=b(5,iti) < b2(1,i-2)=b(2,iti) < b2(2,i-2)=b(4,iti) < b1tilde(1,i-2)=b1(1,i-2) < b1tilde(2,i-2)=-b1(2,i-2) < b2tilde(1,i-2)=b2(1,i-2) < b2tilde(2,i-2)=-b2(2,i-2) < EE(1,2,i-2)=eeold(1,2,iti) < EE(2,1,i-2)=eeold(2,1,iti) < EE(2,2,i-2)=eeold(2,2,iti) < EE(1,1,i-2)=eeold(1,1,iti) < enddo < #endif < #ifdef PARMAT < do i=ivec_start+2,ivec_end+2 < #else < do i=3,nres+1 < #endif 2964,2965c2422,2423 < c write(iout,*) "co jest kurwa", iti, EE(1,1,i),EE(2,1,i), < c & EE(1,2,iti),EE(2,2,i) --- > c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti), > c & EE(1,2,iti),EE(2,2,iti) 3010,3011c2468,2476 < C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) < c write (iout,*) 'mu ',mu(:,i-2),i-2 --- > #ifdef MUOUT > write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1), > & rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2), > & -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2), > & dsqrt(b2(1,i-1)**2+b2(2,i-1)**2) > & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2), > & ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itype2loc(iti)) > #endif > cd write (iout,*) 'mu ',mu(:,i-2) 3298c2763 < cd iti = itortyp(itype(i)) --- > cd iti = itype2loc(itype(i)) 3335d2799 < include 'COMMON.SPLITELE' 3418d2881 < C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition 3420,3421d2882 < if (i.le.1) cycle < C write(iout,*) "tu jest i",i 3423,3434c2884 < C changes suggested by Ana to avoid out of bounds < & .or.((i+4).gt.nres) < & .or.((i-1).le.0) < C end of changes by Ana < & .or. itype(i+2).eq.ntyp1 < & .or. itype(i+3).eq.ntyp1) cycle < if(i.gt.1)then < if(itype(i-1).eq.ntyp1)cycle < end if < if(i.LT.nres-3)then < if (itype(i+4).eq.ntyp1) cycle < end if --- > & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle 3444,3449d2893 < xmedi=mod(xmedi,boxxsize) < if (xmedi.lt.0) xmedi=xmedi+boxxsize < ymedi=mod(ymedi,boxysize) < if (ymedi.lt.0) ymedi=ymedi+boxysize < zmedi=mod(zmedi,boxzsize) < if (zmedi.lt.0) zmedi=zmedi+boxzsize 3456d2899 < if (i.le.1) cycle 3458,3461d2900 < C changes suggested by Ana to avoid out of bounds < & .or.((i+5).gt.nres) < & .or.((i-1).le.0) < C end of changes suggested by Ana 3463,3467c2902 < & .or. itype(i+4).eq.ntyp1 < & .or. itype(i+5).eq.ntyp1 < & .or. itype(i).eq.ntyp1 < & .or. itype(i-1).eq.ntyp1 < & ) cycle --- > & .or. itype(i+4).eq.ntyp1) cycle 3477,3508d2911 < C Return atom into box, boxxsize is size of box in x dimension < c 194 continue < c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize < c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize < C Condition for being inside the proper box < c if ((xmedi.gt.((0.5d0)*boxxsize)).or. < c & (xmedi.lt.((-0.5d0)*boxxsize))) then < c go to 194 < c endif < c 195 continue < c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize < c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize < C Condition for being inside the proper box < c if ((ymedi.gt.((0.5d0)*boxysize)).or. < c & (ymedi.lt.((-0.5d0)*boxysize))) then < c go to 195 < c endif < c 196 continue < c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize < c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize < C Condition for being inside the proper box < c if ((zmedi.gt.((0.5d0)*boxzsize)).or. < c & (zmedi.lt.((-0.5d0)*boxzsize))) then < c go to 196 < c endif < xmedi=mod(xmedi,boxxsize) < if (xmedi.lt.0) xmedi=xmedi+boxxsize < ymedi=mod(ymedi,boxysize) < if (ymedi.lt.0) ymedi=ymedi+boxysize < zmedi=mod(zmedi,boxzsize) < if (zmedi.lt.0) zmedi=zmedi+boxzsize < 3516,3519d2918 < C Loop over all neighbouring boxes < C do xshift=-1,1 < C do yshift=-1,1 < C do zshift=-1,1 3523d2921 < CTU KURWA 3525,3534c2923,2924 < C do i=75,75 < if (i.le.1) cycle < if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 < C changes suggested by Ana to avoid out of bounds < & .or.((i+2).gt.nres) < & .or.((i-1).le.0) < C end of changes by Ana < & .or. itype(i+2).eq.ntyp1 < & .or. itype(i-1).eq.ntyp1 < & ) cycle --- > c do i=7,7 > if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle 3544,3579d2933 < xmedi=mod(xmedi,boxxsize) < if (xmedi.lt.0) xmedi=xmedi+boxxsize < ymedi=mod(ymedi,boxysize) < if (ymedi.lt.0) ymedi=ymedi+boxysize < zmedi=mod(zmedi,boxzsize) < if (zmedi.lt.0) zmedi=zmedi+boxzsize < C xmedi=xmedi+xshift*boxxsize < C ymedi=ymedi+yshift*boxysize < C zmedi=zmedi+zshift*boxzsize < < C Return tom into box, boxxsize is size of box in x dimension < c 164 continue < c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize < c if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize < C Condition for being inside the proper box < c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or. < c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then < c go to 164 < c endif < c 165 continue < c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize < c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize < C Condition for being inside the proper box < c if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or. < c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then < c go to 165 < c endif < c 166 continue < c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize < c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize < cC Condition for being inside the proper box < c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or. < c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then < c go to 166 < c endif < 3582d2935 < C I TU KURWA 3584,3594c2937,2939 < C do j=16,17 < C write (iout,*) i,j < if (j.le.1) cycle < if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 < C changes suggested by Ana to avoid out of bounds < & .or.((j+2).gt.nres) < & .or.((j-1).le.0) < C end of changes by Ana < & .or.itype(j+2).eq.ntyp1 < & .or.itype(j-1).eq.ntyp1 < &) cycle --- > c do j=13,13 > c write (iout,*) 'tu wchodze',i,j,itype(i),itype(j) > if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle 3599,3602d2943 < C enddo ! zshift < C enddo ! yshift < C enddo ! xshift < 3633,3634d2973 < include 'COMMON.SPLITELE' < include 'COMMON.SHIELD' 3670,3742c3009,3011 < C xj=c(1,j)+0.5D0*dxj-xmedi < C yj=c(2,j)+0.5D0*dyj-ymedi < C zj=c(3,j)+0.5D0*dzj-zmedi < xj=c(1,j)+0.5D0*dxj < yj=c(2,j)+0.5D0*dyj < zj=c(3,j)+0.5D0*dzj < xj=mod(xj,boxxsize) < if (xj.lt.0) xj=xj+boxxsize < yj=mod(yj,boxysize) < if (yj.lt.0) yj=yj+boxysize < zj=mod(zj,boxzsize) < if (zj.lt.0) zj=zj+boxzsize < if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" < dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 < xj_safe=xj < yj_safe=yj < zj_safe=zj < isubchap=0 < do xshift=-1,1 < do yshift=-1,1 < do zshift=-1,1 < xj=xj_safe+xshift*boxxsize < yj=yj_safe+yshift*boxysize < zj=zj_safe+zshift*boxzsize < dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 < if(dist_temp.lt.dist_init) then < dist_init=dist_temp < xj_temp=xj < yj_temp=yj < zj_temp=zj < isubchap=1 < endif < enddo < enddo < enddo < if (isubchap.eq.1) then < xj=xj_temp-xmedi < yj=yj_temp-ymedi < zj=zj_temp-zmedi < else < xj=xj_safe-xmedi < yj=yj_safe-ymedi < zj=zj_safe-zmedi < endif < C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC < c 174 continue < c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize < c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize < C Condition for being inside the proper box < c if ((xj.gt.((0.5d0)*boxxsize)).or. < c & (xj.lt.((-0.5d0)*boxxsize))) then < c go to 174 < c endif < c 175 continue < c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize < c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize < C Condition for being inside the proper box < c if ((yj.gt.((0.5d0)*boxysize)).or. < c & (yj.lt.((-0.5d0)*boxysize))) then < c go to 175 < c endif < c 176 continue < c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize < c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize < C Condition for being inside the proper box < c if ((zj.gt.((0.5d0)*boxzsize)).or. < c & (zj.lt.((-0.5d0)*boxzsize))) then < c go to 176 < c endif < C endif !endPBC condintion < C xj=xj-xmedi < C yj=yj-ymedi < C zj=zj-zmedi --- > xj=c(1,j)+0.5D0*dxj-xmedi > yj=c(2,j)+0.5D0*dyj-ymedi > zj=c(3,j)+0.5D0*dzj-zmedi 3744,3747d3012 < < sss=sscale(sqrt(rij)) < sssgrad=sscagrad(sqrt(rij)) < c if (sss.gt.0.0d0) then 3763c3028 < evdwij=(ev1+ev2) --- > evdwij=ev1+ev2 3766,3767c3031 < C MARYSIA < C eesij=(el1+el2) --- > eesij=el1+el2 3770,3780d3033 < if (shield_mode.gt.0) then < C fac_shield(i)=0.4 < C fac_shield(j)=0.6 < el1=el1*fac_shield(i)**2*fac_shield(j)**2 < el2=el2*fac_shield(i)**2*fac_shield(j)**2 < eesij=(el1+el2) < ees=ees+eesij < else < fac_shield(i)=1.0 < fac_shield(j)=1.0 < eesij=(el1+el2) 3782,3783c3035 < endif < evdw1=evdw1+evdwij*sss --- > evdw1=evdw1+evdwij 3793,3794c3045 < write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, < &fac_shield(i),fac_shield(j) --- > write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij 3801c3052 < facvdw=-6*rrmij*(ev1+evdwij)*sss --- > facvdw=-6*rrmij*(ev1+evdwij) 3807d3057 < 3814,3882d3063 < if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. < & (shield_mode.gt.0)) then < C print *,i,j < do ilist=1,ishield_list(i) < iresshield=shield_list(ilist,i) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i) < & *2.0 < gshieldx(k,iresshield)=gshieldx(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 < gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield < C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) < C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) < C if (iresshield.gt.i) then < C do ishi=i+1,iresshield-1 < C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield < C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) < C < C enddo < C else < C do ishi=iresshield,i < C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield < C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) < C < C enddo < C endif < enddo < enddo < do ilist=1,ishield_list(j) < iresshield=shield_list(ilist,j) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) < & *2.0 < gshieldx(k,iresshield)=gshieldx(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 < gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield < < C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) < C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) < C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) < C if (iresshield.gt.j) then < C do ishi=j+1,iresshield-1 < C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield < C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) < C < C enddo < C else < C do ishi=iresshield,j < C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield < C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) < C enddo < C endif < enddo < enddo < < do k=1,3 < gshieldc(k,i)=gshieldc(k,i)+ < & grad_shield(k,i)*eesij/fac_shield(i)*2.0 < gshieldc(k,j)=gshieldc(k,j)+ < & grad_shield(k,j)*eesij/fac_shield(j)*2.0 < gshieldc(k,i-1)=gshieldc(k,i-1)+ < & grad_shield(k,i)*eesij/fac_shield(i)*2.0 < gshieldc(k,j-1)=gshieldc(k,j-1)+ < & grad_shield(k,j)*eesij/fac_shield(j)*2.0 < < enddo < endif 3889d3069 < C print *,"before", gelc_long(1,i), gelc_long(1,j) 3892d3071 < C & +grad_shield(k,j)*eesij/fac_shield(j) 3894,3898d3072 < C & +grad_shield(k,i)*eesij/fac_shield(i) < C gelc_long(k,i-1)=gelc_long(k,i-1) < C & +grad_shield(k,i)*eesij/fac_shield(i) < C gelc_long(k,j-1)=gelc_long(k,j-1) < C & +grad_shield(k,j)*eesij/fac_shield(j) 3900,3901d3073 < C print *,"bafter", gelc_long(1,i), gelc_long(1,j) < 3910,3918c3082,3084 < if (sss.gt.0.0) then < ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj < ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj < ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj < else < ggg(1)=0.0 < ggg(2)=0.0 < ggg(3)=0.0 < endif --- > ggg(1)=facvdw*xj > ggg(2)=facvdw*yj > ggg(3)=facvdw*zj 3938,3940c3104,3105 < C MARYSIA < facvdw=(ev1+evdwij)*sss < facel=(el1+eesij) --- > facvdw=ev1+evdwij > facel=el1+eesij 3950d3114 < C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j) 3952d3115 < C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j) 3954d3116 < C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j) 3974,3976c3136,3138 < ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj < ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj < ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj --- > ggg(1)=facvdw*xj > ggg(2)=facvdw*yj > ggg(3)=facvdw*zj 3997,3998c3159 < ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* < & fac_shield(i)**2*fac_shield(j)**2 --- > ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) 4014d3174 < C print *,"before22", gelc_long(1,i), gelc_long(1,j) 4017,4019c3177,3178 < & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) < & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) < & *fac_shield(i)**2*fac_shield(j)**2 --- > & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) > & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) 4021,4023c3180,3181 < & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) < & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) < & *fac_shield(i)**2*fac_shield(j)**2 --- > & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) > & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) 4027,4030d3184 < C print *,"before33", gelc_long(1,i), gelc_long(1,j) < < C MARYSIA < c endif !sscale 4055d3208 < c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l 4058c3211 < c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) --- > c write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i 4061c3214 < c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) --- > c write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j 4083c3236 < cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33 --- > cd & i,itype2loc(itype(i)),j,itype2loc(itype(j)),a22,a23,a32,a33 4225a3379 > c if ((i.eq.8).and.(j.eq.14)) then 4229,4286d3382 < if (shield_mode.eq.0) then < fac_shield(i)=1.0 < fac_shield(j)=1.0 < C else < C fac_shield(i)=0.4 < C fac_shield(j)=0.6 < endif < eel_loc_ij=eel_loc_ij < & *fac_shield(i)*fac_shield(j) < C Now derivative over eel_loc < if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. < & (shield_mode.gt.0)) then < C print *,i,j < < do ilist=1,ishield_list(i) < iresshield=shield_list(ilist,i) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij < & /fac_shield(i) < C & *2.0 < gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) < gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) < & +rlocshield < enddo < enddo < do ilist=1,ishield_list(j) < iresshield=shield_list(ilist,j) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij < & /fac_shield(j) < C & *2.0 < gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) < gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1) < & +rlocshield < < enddo < enddo < < do k=1,3 < gshieldc_ll(k,i)=gshieldc_ll(k,i)+ < & grad_shield(k,i)*eel_loc_ij/fac_shield(i) < gshieldc_ll(k,j)=gshieldc_ll(k,j)+ < & grad_shield(k,j)*eel_loc_ij/fac_shield(j) < gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ < & grad_shield(k,i)*eel_loc_ij/fac_shield(i) < gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ < & grad_shield(k,j)*eel_loc_ij/fac_shield(j) < enddo < endif < < < c write (iout,*) 'i',i,' j',j,itype(i),itype(j), < c & ' eel_loc_ij',eel_loc_ij < C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) 4289c3385 < geel_loc_ij=(a22*gmuij1(1) --- > geel_loc_ij=a22*gmuij1(1) 4292,4293c3388 < & +a33*gmuij1(4)) < & *fac_shield(i)*fac_shield(j) --- > & +a33*gmuij1(4) 4302,4305c3397 < geel_loc_ij= < & a22*gmuij2(1) < & +a23*gmuij2(2) < & +a32*gmuij2(3) --- > geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3) 4309,4314c3401 < & *fac_shield(i)*fac_shield(j) < < c Derivative over j residue < geel_loc_ji=a22*gmuji1(1) < & +a23*gmuji1(2) < & +a32*gmuji1(3) --- > geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3) 4320c3407 < gloc(nphi+j,icg)=gloc(nphi+j,icg)+ --- > gloc(nphi+j,icg)=gloc(nphi+j,icg)+ 4322,4327c3409 < & *fac_shield(i)*fac_shield(j) < < geel_loc_ji= < & +a22*gmuji2(1) < & +a23*gmuji2(2) < & +a32*gmuji2(3) --- > geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3) 4334d3415 < & *fac_shield(i)*fac_shield(j) 4340,4342c3421 < c if (eel_loc_ij.ne.0) < c & write (iout,'(a4,2i4,8f9.5)')'chuj', < c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) --- > c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) 4348,4351c3427,3428 < & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) < & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) < & *fac_shield(i)*fac_shield(j) < --- > & a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) > & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j) 4353,4355c3430,3431 < & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) < & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) < & *fac_shield(i)*fac_shield(j) --- > & a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) > & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j) 4358,4360c3434,3435 < ggg(l)=(agg(l,1)*muij(1)+ < & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) < & *fac_shield(i)*fac_shield(j) --- > ggg(l)=agg(l,1)*muij(1)+ > & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4) 4374,4389c3449,3456 < gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ < & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) < & *fac_shield(i)*fac_shield(j) < < gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ < & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) < & *fac_shield(i)*fac_shield(j) < < gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ < & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) < & *fac_shield(i)*fac_shield(j) < < gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ < & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) < & *fac_shield(i)*fac_shield(j) < --- > gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ > & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4) > gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+ > & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4) > gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+ > & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4) > gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+ > & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4) 4390a3458 > c endif 4468,4475d3535 < if (shield_mode.eq.0) then < fac_shield(i)=1.0d0 < fac_shield(j)=1.0d0 < else < ees0plist(num_conti,i)=j < C fac_shield(i)=0.4d0 < C fac_shield(j)=0.6d0 < endif 4477d3536 < & *fac_shield(i)*fac_shield(j) 4479d3537 < & *fac_shield(i)*fac_shield(j) 4547,4548d3604 < & *fac_shield(i)*fac_shield(j) < 4552,4553d3607 < & *fac_shield(i)*fac_shield(j) < 4555,4556d3608 < & *fac_shield(i)*fac_shield(j) < 4560,4561d3611 < & *fac_shield(i)*fac_shield(j) < 4565,4566d3614 < & *fac_shield(i)*fac_shield(j) < 4568,4569d3615 < & *fac_shield(i)*fac_shield(j) < 4621d3666 < include 'COMMON.SHIELD' 4662,4668d3706 < if (shield_mode.eq.0) then < fac_shield(i)=1.0 < fac_shield(j)=1.0 < C else < C fac_shield(i)=0.4 < C fac_shield(j)=0.6 < endif 4670,4672d3707 < & *fac_shield(i)*fac_shield(j) < eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) < & *fac_shield(i)*fac_shield(j) 4676d3710 < & *fac_shield(i)*fac_shield(j) 4679,4709d3712 < & *fac_shield(i)*fac_shield(j) < < < C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') < C Derivatives in shield mode < if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. < & (shield_mode.gt.0)) then < C print *,i,j < < do ilist=1,ishield_list(i) < iresshield=shield_list(ilist,i) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i) < C & *2.0 < gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i) < gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) < & +rlocshield < enddo < enddo < do ilist=1,ishield_list(j) < iresshield=shield_list(ilist,j) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j) < C & *2.0 < gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j) < gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1) < & +rlocshield 4711,4726c3714,3715 < enddo < enddo < < do k=1,3 < gshieldc_t3(k,i)=gshieldc_t3(k,i)+ < & grad_shield(k,i)*eello_t3/fac_shield(i) < gshieldc_t3(k,j)=gshieldc_t3(k,j)+ < & grad_shield(k,j)*eello_t3/fac_shield(j) < gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+ < & grad_shield(k,i)*eello_t3/fac_shield(i) < gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+ < & grad_shield(k,j)*eello_t3/fac_shield(j) < enddo < endif < < C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') --- > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) 4735d3723 < & *fac_shield(i)*fac_shield(j) 4742d3729 < & *fac_shield(i)*fac_shield(j) 4756,4757d3742 < & *fac_shield(i)*fac_shield(j) < 4765d3749 < & *fac_shield(i)*fac_shield(j) 4773d3756 < & *fac_shield(i)*fac_shield(j) 4781d3763 < & *fac_shield(i)*fac_shield(j) 4802d3783 < include 'COMMON.SHIELD' 4903,4961d3883 < if (shield_mode.eq.0) then < fac_shield(i)=1.0 < fac_shield(j)=1.0 < C else < C fac_shield(i)=0.6 < C fac_shield(j)=0.4 < endif < eello_turn4=eello_turn4-(s1+s2+s3) < & *fac_shield(i)*fac_shield(j) < eello_t4=-(s1+s2+s3) < & *fac_shield(i)*fac_shield(j) < c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) < if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)') < & 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3 < C Now derivative over shield: < if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. < & (shield_mode.gt.0)) then < C print *,i,j < < do ilist=1,ishield_list(i) < iresshield=shield_list(ilist,i) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i) < C & *2.0 < gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i) < gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) < & +rlocshield < enddo < enddo < do ilist=1,ishield_list(j) < iresshield=shield_list(ilist,j) < do k=1,3 < rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j) < C & *2.0 < gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ < & rlocshield < & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j) < gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) < & +rlocshield < < enddo < enddo < < do k=1,3 < gshieldc_t4(k,i)=gshieldc_t4(k,i)+ < & grad_shield(k,i)*eello_t4/fac_shield(i) < gshieldc_t4(k,j)=gshieldc_t4(k,j)+ < & grad_shield(k,j)*eello_t4/fac_shield(j) < gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+ < & grad_shield(k,i)*eello_t4/fac_shield(i) < gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+ < & grad_shield(k,j)*eello_t4/fac_shield(j) < enddo < endif < < < 4963,4966c3885 < < < cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), < cd & ' eello_turn4_num',8*eello_turn4_num --- > eello_turn4=eello_turn4-(s1+s2+s3) 4970d3888 < & *fac_shield(i)*fac_shield(j) 4973,4974d3890 < & *fac_shield(i)*fac_shield(j) < 4977,4978d3892 < & *fac_shield(i)*fac_shield(j) < 4994d3907 < & *fac_shield(i)*fac_shield(j) 5003d3915 < & *fac_shield(i)*fac_shield(j) 5015d3926 < & *fac_shield(i)*fac_shield(j) 5035d3945 < & *fac_shield(i)*fac_shield(j) 5054d3963 < & *fac_shield(i)*fac_shield(j) 5069d3977 < & *fac_shield(i)*fac_shield(j) 5084d3991 < & *fac_shield(i)*fac_shield(j) 5100d4006 < & *fac_shield(i)*fac_shield(j) 5161,5163d4066 < C do xshift=-1,1 < C do yshift=-1,1 < C do zshift=-1,1 5170,5203c4073 < C Return atom into box, boxxsize is size of box in x dimension < c 134 continue < c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize < c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize < C Condition for being inside the proper box < c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. < c & (xi.lt.((xshift-0.5d0)*boxxsize))) then < c go to 134 < c endif < c 135 continue < c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize < c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize < C Condition for being inside the proper box < c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. < c & (yi.lt.((yshift-0.5d0)*boxysize))) then < c go to 135 < c c endif < c 136 continue < c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize < c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize < cC Condition for being inside the proper box < c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. < c & (zi.lt.((zshift-0.5d0)*boxzsize))) then < c go to 136 < c endif < xi=mod(xi,boxxsize) < if (xi.lt.0) xi=xi+boxxsize < yi=mod(yi,boxysize) < if (yi.lt.0) yi=yi+boxysize < zi=mod(zi,boxzsize) < if (zi.lt.0) zi=zi+boxzsize < C xi=xi+xshift*boxxsize < C yi=yi+yshift*boxysize < C zi=zi+zshift*boxzsize --- > 5214,5280c4084,4086 < xj=c(1,j) < yj=c(2,j) < zj=c(3,j) < c 174 continue < c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize < c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize < C Condition for being inside the proper box < c if ((xj.gt.((0.5d0)*boxxsize)).or. < c & (xj.lt.((-0.5d0)*boxxsize))) then < c go to 174 < c endif < c 175 continue < c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize < c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize < cC Condition for being inside the proper box < c if ((yj.gt.((0.5d0)*boxysize)).or. < c & (yj.lt.((-0.5d0)*boxysize))) then < c go to 175 < c endif < c 176 continue < c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize < c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize < C Condition for being inside the proper box < c if ((zj.gt.((0.5d0)*boxzsize)).or. < c & (zj.lt.((-0.5d0)*boxzsize))) then < c go to 176 < xj=mod(xj,boxxsize) < if (xj.lt.0) xj=xj+boxxsize < yj=mod(yj,boxysize) < if (yj.lt.0) yj=yj+boxysize < zj=mod(zj,boxzsize) < if (zj.lt.0) zj=zj+boxzsize < dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < xj_safe=xj < yj_safe=yj < zj_safe=zj < subchap=0 < do xshift=-1,1 < do yshift=-1,1 < do zshift=-1,1 < xj=xj_safe+xshift*boxxsize < yj=yj_safe+yshift*boxysize < zj=zj_safe+zshift*boxzsize < dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < if(dist_temp.lt.dist_init) then < dist_init=dist_temp < xj_temp=xj < yj_temp=yj < zj_temp=zj < subchap=1 < endif < enddo < enddo < enddo < if (subchap.eq.1) then < xj=xj_temp-xi < yj=yj_temp-yi < zj=zj_temp-zi < else < xj=xj_safe-xi < yj=yj_safe-yi < zj=zj_safe-zi < endif < c c endif < C xj=xj-xi < C yj=yj-yi < C zj=zj-zi --- > xj=c(1,j)-xi > yj=c(2,j)-yi > zj=c(3,j)-zi 5282d4087 < 5333,5335d4137 < C enddo !zshift < C enddo !yshift < C enddo !xshift 5356d4157 < include 'COMMON.SPLITELE' 5360d4160 < c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' 5363,5365d4162 < C do xshift=-1,1 < C do yshift=-1,1 < C do zshift=-1,1 5372,5392d4168 < xi=mod(xi,boxxsize) < if (xi.lt.0) xi=xi+boxxsize < yi=mod(yi,boxysize) < if (yi.lt.0) yi=yi+boxysize < zi=mod(zi,boxzsize) < if (zi.lt.0) zi=zi+boxzsize < c xi=xi+xshift*boxxsize < c yi=yi+yshift*boxysize < c zi=zi+zshift*boxzsize < c print *,xi,yi,zi,'polozenie i' < C Return atom into box, boxxsize is size of box in x dimension < c 134 continue < c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize < c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize < C Condition for being inside the proper box < c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. < c & (xi.lt.((xshift-0.5d0)*boxxsize))) then < c go to 134 < c endif < c 135 continue < c print *,xi,boxxsize,"pierwszy" 5394,5408d4169 < c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize < c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize < C Condition for being inside the proper box < c if ((yi.gt.((yshift+0.5d0)*boxysize)).or. < c & (yi.lt.((yshift-0.5d0)*boxysize))) then < c go to 135 < c endif < c 136 continue < c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize < c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize < C Condition for being inside the proper box < c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. < c & (zi.lt.((zshift-0.5d0)*boxzsize))) then < c go to 136 < c endif 5419,5484c4180,4182 < xj=c(1,j) < yj=c(2,j) < zj=c(3,j) < xj=mod(xj,boxxsize) < if (xj.lt.0) xj=xj+boxxsize < yj=mod(yj,boxysize) < if (yj.lt.0) yj=yj+boxysize < zj=mod(zj,boxzsize) < if (zj.lt.0) zj=zj+boxzsize < c 174 continue < c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize < c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize < C Condition for being inside the proper box < c if ((xj.gt.((0.5d0)*boxxsize)).or. < c & (xj.lt.((-0.5d0)*boxxsize))) then < c go to 174 < c endif < c 175 continue < c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize < c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize < cC Condition for being inside the proper box < c if ((yj.gt.((0.5d0)*boxysize)).or. < c & (yj.lt.((-0.5d0)*boxysize))) then < c go to 175 < c endif < c 176 continue < c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize < c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize < C Condition for being inside the proper box < c if ((zj.gt.((0.5d0)*boxzsize)).or. < c & (zj.lt.((-0.5d0)*boxzsize))) then < c go to 176 < c endif < CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE < dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < xj_safe=xj < yj_safe=yj < zj_safe=zj < subchap=0 < do xshift=-1,1 < do yshift=-1,1 < do zshift=-1,1 < xj=xj_safe+xshift*boxxsize < yj=yj_safe+yshift*boxysize < zj=zj_safe+zshift*boxzsize < dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 < if(dist_temp.lt.dist_init) then < dist_init=dist_temp < xj_temp=xj < yj_temp=yj < zj_temp=zj < subchap=1 < endif < enddo < enddo < enddo < if (subchap.eq.1) then < xj=xj_temp-xi < yj=yj_temp-yi < zj=zj_temp-zi < else < xj=xj_safe-xi < yj=yj_safe-yi < zj=zj_safe-zi < endif < c print *,xj,yj,zj,'polozenie j' --- > xj=c(1,j)-xi > yj=c(2,j)-yi > zj=c(3,j)-zi 5486,5491d4183 < c print *,rrij < sss=sscale(1.0d0/(dsqrt(rrij))) < c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz' < c if (sss.eq.0) print *,'czasem jest OK' < if (sss.le.0.0d0) cycle < sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) 5498c4190 < evdw2_14=evdw2_14+(e1+e2)*sss --- > evdw2_14=evdw2_14+e1+e2 5501c4193 < evdw2=evdw2+evdwij*sss --- > evdw2=evdw2+evdwij 5508,5509c4200 < fac=-(evdwij+e1)*rrij*sss < fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon --- > fac=-(evdwij+e1)*rrij 5544,5545c4235 < c endif !endif for sscale cutoff < enddo ! j --- > enddo 5549,5551d4238 < c enddo !zshift < c enddo !yshift < c enddo !xshift 5583d4269 < include 'COMMON.CONTROL' 5586,5589d4271 < do i=1,3 < ggg(i)=0.0d0 < enddo < C write (iout,*) ,"link_end",link_end,constr_dist 5606,5607c4288 < c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, < c & dhpb(i),dhpb1(i),forcon(i) --- > cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj 5610,5616c4291 < C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. < C & iabs(itype(jjj)).eq.1) then < cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then < C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds < if (.not.dyn_ss .and. i.le.nss) then < C 15/02/13 CC dynamic SSbond - additional check < if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. --- > if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. 5620d4294 < endif 5622,5663d4295 < cd & ' waga=',waga,' fac=',fac < else if (ii.gt.nres .and. jj.gt.nres) then < c Restraints from contact prediction < dd=dist(ii,jj) < if (constr_dist.eq.11) then < ehpb=ehpb+fordepth(i)**4.0d0 < & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) < fac=fordepth(i)**4.0d0 < & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd < if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, < & ehpb,fordepth(i),dd < else < if (dhpb1(i).gt.0.0d0) then < ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) < fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd < c write (iout,*) "beta nmr", < c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) < else < dd=dist(ii,jj) < rdis=dd-dhpb(i) < C Get the force constant corresponding to this distance. < waga=forcon(i) < C Calculate the contribution to energy. < ehpb=ehpb+waga*rdis*rdis < c write (iout,*) "beta reg",dd,waga*rdis*rdis < C < C Evaluate gradient. < C < fac=waga*rdis/dd < endif < endif < do j=1,3 < ggg(j)=fac*(c(j,jj)-c(j,ii)) < enddo < do j=1,3 < ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) < ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) < enddo < do k=1,3 < ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) < ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) < enddo 5667,5682c4299,4300 < dd=dist(ii,jj) < if (constr_dist.eq.11) then < ehpb=ehpb+fordepth(i)**4.0d0 < & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) < fac=fordepth(i)**4.0d0 < & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd < if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, < & ehpb,fordepth(i),dd < else < if (dhpb1(i).gt.0.0d0) then < ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) < fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd < c write (iout,*) "alph nmr", < c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) < else < rdis=dd-dhpb(i) --- > dd=dist(ii,jj) > rdis=dd-dhpb(i) 5684c4302 < waga=forcon(i) --- > waga=forcon(i) 5686,5687c4304 < ehpb=ehpb+waga*rdis*rdis < c write (iout,*) "alpha reg",dd,waga*rdis*rdis --- > ehpb=ehpb+waga*rdis*rdis 5691,5696c4308,4313 < fac=waga*rdis/dd < endif < endif < do j=1,3 < ggg(j)=fac*(c(j,jj)-c(j,ii)) < enddo --- > fac=waga*rdis/dd > cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, > cd & ' waga=',waga,' fac=',fac > do j=1,3 > ggg(j)=fac*(c(j,jj)-c(j,ii)) > enddo 5700c4317 < if (iii.lt.ii) then --- > if (iii.lt.ii) then 5705c4322 < endif --- > endif 5711,5714c4328,4331 < do k=1,3 < ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) < ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) < enddo --- > do k=1,3 > ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) > ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) > enddo 5717c4334 < if (constr_dist.ne.11) ehpb=0.5D0*ehpb --- > ehpb=0.5D0*ehpb 5830,5844c4447,4455 < if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle < c estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) < c do j=1,3 < c gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) < c & *dc(j,i-1)/vbld(i) < c enddo < c if (energy_dec) write(iout,*) < c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) < c else < C Checking if it involves dummy (NH3+ or COO-) group < if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then < C YES vbldpDUM is the equlibrium length of spring for Dummy atom < diff = vbld(i)-vbldpDUM < else < C NO vbldp0 is the equlibrium lenght of spring for peptide group --- > if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then > estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) > do j=1,3 > gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) > & *dc(j,i-1)/vbld(i) > enddo > if (energy_dec) write(iout,*) > & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) > else 5846,5847c4457 < endif < if (energy_dec) write (iout,'(a7,i5,4f7.3)') --- > if (energy_dec) write (iout,*) 5854c4464 < c endif --- > endif 5866c4476 < if (energy_dec) write (iout,*) --- > if (energy_dec) write (iout,*) 5908c4518 < subroutine ebend(etheta,ethetacnstr) --- > subroutine ebend(etheta) 5925d4534 < include 'COMMON.TORCNSTR' 5936,5937c4545 < if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 < & .or.itype(i).eq.ntyp1) cycle --- > if (itype(i-1).eq.ntyp1) cycle 5954c4562 < if (i.gt.3 .and. itype(i-3).ne.ntyp1) then --- > if (i.gt.3 .and. itype(i-2).ne.ntyp1) then 5967c4575 < if (i.lt.nres .and. itype(i+1).ne.ntyp1) then --- > if (i.lt.nres .and. itype(i).ne.ntyp1) then 5975d4582 < #endif 5976a4584 > #endif 5994d4601 < c write(iout,*) 'chuj tu', y(k),z(k) 6031,6032c4638,4639 < if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)') < & 'ebend',i,ethetai,theta(i),itype(i) --- > if (energy_dec) write (iout,'(a6,i5,0pf7.3)') > & 'ebend',i,ethetai 6037,6064d4643 < ethetacnstr=0.0d0 < C print *,ithetaconstr_start,ithetaconstr_end,"TU" < do i=ithetaconstr_start,ithetaconstr_end < itheta=itheta_constr(i) < thetiii=theta(itheta) < difi=pinorm(thetiii-theta_constr0(i)) < if (difi.gt.theta_drange(i)) then < difi=difi-theta_drange(i) < ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 < gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) < & +for_thet_constr(i)*difi**3 < else if (difi.lt.-drange(i)) then < difi=difi+drange(i) < ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 < gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) < & +for_thet_constr(i)*difi**3 < else < difi=0.0 < endif < if (energy_dec) then < write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", < & i,itheta,rad2deg*thetiii, < & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), < & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, < & gloc(itheta+nphi-2,icg) < endif < enddo < 6081,6082c4660 < C the distributioni. < ccc write (iout,*) thetai,thet_pred_mean --- > C the distribution. 6112d4689 < C write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0 6140d4716 < C write (iout,*) 'termexp',termexp,termm,termpre,i 6181c4757 < subroutine ebend(etheta,ethetacnstr) --- > subroutine ebend(etheta) 6200d4775 < include 'COMMON.TORCNSTR' 6208,6211c4783 < c print *,i,itype(i-1),itype(i),itype(i-2) < if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 < & .or.itype(i).eq.ntyp1) cycle < C print *,i,theta(i) --- > if (itype(i-1).eq.ntyp1) cycle 6223,6224c4795 < C print *,ethetai < if (i.gt.3 .and. itype(i-3).ne.ntyp1) then --- > if (i.gt.3 .and. itype(i-2).ne.ntyp1) then 6238a4810 > ityp1=nthetyp+1 6240d4811 < ityp1=ithetyp((itype(i-2))) 6245c4816 < if (i.lt.nres .and. itype(i+1).ne.ntyp1) then --- > if (i.lt.nres .and. itype(i).ne.ntyp1) then 6260c4831 < ityp3=ithetyp((itype(i))) --- > ityp3=nthetyp+1 6310d4880 < C print *,ethetai 6331d4900 < C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) 6334,6337d4902 < C print *,"cosph1", (cosph1(k), k=1,nsingle) < C print *,"cosph2", (cosph2(k), k=1,nsingle) < C print *,"sinph1", (sinph1(k), k=1,nsingle) < C print *,"sinph2", (sinph2(k), k=1,nsingle) 6340d4904 < C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) 6376d4939 < C print *,ethetai 6385c4948 < gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai --- > gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg) 6387,6415d4949 < C now constrains < ethetacnstr=0.0d0 < C print *,ithetaconstr_start,ithetaconstr_end,"TU" < do i=ithetaconstr_start,ithetaconstr_end < itheta=itheta_constr(i) < thetiii=theta(itheta) < difi=pinorm(thetiii-theta_constr0(i)) < if (difi.gt.theta_drange(i)) then < difi=difi-theta_drange(i) < ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 < gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) < & +for_thet_constr(i)*difi**3 < else if (difi.lt.-drange(i)) then < difi=difi+drange(i) < ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 < gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) < & +for_thet_constr(i)*difi**3 < else < difi=0.0 < endif < if (energy_dec) then < write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", < & i,itheta,rad2deg*thetiii, < & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), < & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, < & gloc(itheta+nphi-2,icg) < endif < enddo < 7180,7182c5714,5716 < & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle < itori=itortyp(itype(i-2)) < itori1=itortyp(itype(i-1)) --- > & .or. itype(i).eq.ntyp1) cycle > itori=itortyp(itype(i-2)) > itori1=itortyp(itype(i-1)) 7235,7236c5769,5770 < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 --- > edihcnstr=edihcnstr+0.25d0*ftors*difi**4 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 7239,7240c5773,5774 < edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 --- > edihcnstr=edihcnstr+0.25d0*ftors*difi**4 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 7276,7284c5810,5811 < C ANY TWO ARE DUMMY ATOMS in row CYCLE < c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. < c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. < c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle < if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 < & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle < C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF < C For introducing the NH3+ and COO- group please check the etor_d for reference < C and guidance --- > if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 > & .or. itype(i).eq.ntyp1) cycle 7346,7347c5873,5874 < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 --- > edihcnstr=edihcnstr+0.25d0*ftors*difi**4 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 7350,7351c5877,5878 < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 --- > edihcnstr=edihcnstr+0.25d0*ftors*difi**4 > gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 7355,7360c5882,5884 < if (energy_dec) then < write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc", < & i,itori,rad2deg*phii, < & rad2deg*phi0(i), rad2deg*drange(i), < & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) < endif --- > cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, > cd & rad2deg*phi0(i), rad2deg*drange(i), > cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) 7388,7396c5912,5913 < C ANY TWO ARE DUMMY ATOMS in row CYCLE < C if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. < C & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or. < C & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1)) .or. < C & ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle < if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. < & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or. < & (itype(i+1).eq.ntyp1)) cycle < C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF --- > if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 > & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle 7406,7414d5922 < C Iblock=2 Proline type < C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT < C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO- < C if (itype(i+1).eq.ntyp1) iblock=3 < C The problem of NH3+ group can be resolved by adding new parameters please note if there < C IS or IS NOT need for this < C IF Yes uncomment below and add to parmread.F appropriate changes and to v1cij and so on < C is (itype(i-3).eq.ntyp1) ntblock=2 < C ntblock is N-terminal blocking group 7418,7420d5925 < C Example of changes for NH3+ blocking group < C do j=1,ntermd_1(itori,itori1,itori2,iblock,ntblock) < C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock) 7458,7703d5962 < C---------------------------------------------------------------------------------- < C The rigorous attempt to derive energy function < subroutine etor_kcc(etors,edihcnstr) < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.VAR' < include 'COMMON.GEO' < include 'COMMON.LOCAL' < include 'COMMON.TORSION' < include 'COMMON.INTERACT' < include 'COMMON.DERIV' < include 'COMMON.CHAIN' < include 'COMMON.NAMES' < include 'COMMON.IOUNITS' < include 'COMMON.FFIELD' < include 'COMMON.TORCNSTR' < include 'COMMON.CONTROL' < logical lprn < c double precision thybt1(maxtermkcc),thybt2(maxtermkcc) < C Set lprn=.true. for debugging < lprn=.false. < c lprn=.true. < C print *,"wchodze kcc" < if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode < if (tor_mode.ne.2) then < etors=0.0D0 < endif < do i=iphi_start,iphi_end < C ANY TWO ARE DUMMY ATOMS in row CYCLE < c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. < c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. < c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle < if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 < & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle < itori=itortyp_kcc(itype(i-2)) < itori1=itortyp_kcc(itype(i-1)) < phii=phi(i) < glocig=0.0D0 < glocit1=0.0d0 < glocit2=0.0d0 < sumnonchebyshev=0.0d0 < sumchebyshev=0.0d0 < C to avoid multiple devision by 2 < c theti22=0.5d0*theta(i) < C theta 12 is the theta_1 /2 < C theta 22 is theta_2 /2 < c theti12=0.5d0*theta(i-1) < C and appropriate sinus function < sinthet1=dsin(theta(i-1)) < sinthet2=dsin(theta(i)) < costhet1=dcos(theta(i-1)) < costhet2=dcos(theta(i)) < c Cosines of halves thetas < costheti12=0.5d0*(1.0d0+costhet1) < costheti22=0.5d0*(1.0d0+costhet2) < C to speed up lets store its mutliplication < sint1t2=sinthet2*sinthet1 < sint1t2n=1.0d0 < C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma) < C +d_n*sin(n*gamma)) * < C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2))) < C we have two sum 1) Non-Chebyshev which is with n and gamma < etori=0.0d0 < do j=1,nterm_kcc(itori,itori1) < < nval=nterm_kcc_Tb(itori,itori1) < v1ij=v1_kcc(j,itori,itori1) < v2ij=v2_kcc(j,itori,itori1) < c write (iout,*) "i",i," j",j," v1",v1ij," v2",v2ij < C v1ij is c_n and d_n in euation above < cosphi=dcos(j*phii) < sinphi=dsin(j*phii) < sint1t2n1=sint1t2n < sint1t2n=sint1t2n*sint1t2 < sumth1tyb1=tschebyshev(1,nval,v11_chyb(1,j,itori,itori1), < & costheti12) < gradth1tyb1=-0.5d0*sinthet1*gradtschebyshev(0,nval-1, < & v11_chyb(1,j,itori,itori1),costheti12) < c write (iout,*) "v11",(v11_chyb(k,j,itori,itori1),k=1,nval), < c & " sumth1tyb1",sumth1tyb1," gradth1tyb1",gradth1tyb1 < sumth2tyb1=tschebyshev(1,nval,v21_chyb(1,j,itori,itori1), < & costheti22) < gradth2tyb1=-0.5d0*sinthet2*gradtschebyshev(0,nval-1, < & v21_chyb(1,j,itori,itori1),costheti22) < c write (iout,*) "v21",(v21_chyb(k,j,itori,itori1),k=1,nval), < c & " sumth2tyb1",sumth2tyb1," gradth2tyb1",gradth2tyb1 < sumth1tyb2=tschebyshev(1,nval,v12_chyb(1,j,itori,itori1), < & costheti12) < gradth1tyb2=-0.5d0*sinthet1*gradtschebyshev(0,nval-1, < & v12_chyb(1,j,itori,itori1),costheti12) < c write (iout,*) "v12",(v12_chyb(k,j,itori,itori1),k=1,nval), < c & " sumth1tyb2",sumth1tyb2," gradth1tyb2",gradth1tyb2 < sumth2tyb2=tschebyshev(1,nval,v22_chyb(1,j,itori,itori1), < & costheti22) < gradth2tyb2=-0.5d0*sinthet2*gradtschebyshev(0,nval-1, < & v22_chyb(1,j,itori,itori1),costheti22) < c write (iout,*) "v22",(v22_chyb(k,j,itori,itori1),k=1,nval), < c & " sumth2tyb2",sumth2tyb2," gradth2tyb2",gradth2tyb2 < C etors=etors+sint1t2n*(v1ij*cosphi+v2ij*sinphi) < C if (energy_dec) etors_ii=etors_ii+ < C & v1ij*cosphi+v2ij*sinphi < C glocig is the gradient local i site in gamma < actval1=v1ij*cosphi*(1.0d0+sumth1tyb1+sumth2tyb1) < actval2=v2ij*sinphi*(1.0d0+sumth1tyb2+sumth2tyb2) < etori=etori+sint1t2n*(actval1+actval2) < glocig=glocig+ < & j*sint1t2n*(v2ij*cosphi*(1.0d0+sumth1tyb2+sumth2tyb2) < & -v1ij*sinphi*(1.0d0+sumth1tyb1+sumth2tyb1)) < C now gradient over theta_1 < glocit1=glocit1+ < & j*sint1t2n1*costhet1*sinthet2*(actval1+actval2)+ < & sint1t2n*(v1ij*cosphi*gradth1tyb1+v2ij*sinphi*gradth1tyb2) < glocit2=glocit2+ < & j*sint1t2n1*sinthet1*costhet2*(actval1+actval2)+ < & sint1t2n*(v1ij*cosphi*gradth2tyb1+v2ij*sinphi*gradth2tyb2) < < C now the Czebyshev polinominal sum < c do k=1,nterm_kcc_Tb(itori,itori1) < c thybt1(k)=v1_chyb(k,j,itori,itori1) < c thybt2(k)=v2_chyb(k,j,itori,itori1) < C thybt1(k)=0.0 < C thybt2(k)=0.0 < c enddo < C print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2, < C & gradtschebyshev < C & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1), < C & dcos(theti22)**2), < C & dsin(theti22) < < C now overal sumation < C print *,"sumnon", sumnonchebyshev,sumth1thyb+sumth2thyb < enddo ! j < etors=etors+etori < C derivative over gamma < gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig < C derivative over theta1 < gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1 < C now derivative over theta2 < gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2 < if (lprn) < & write (iout,*) i-2,i-1,itype(i-2),itype(i-1),itori,itori1, < & theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori < enddo < C gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci < ! 6/20/98 - dihedral angle constraints < if (tor_mode.ne.2) then < edihcnstr=0.0d0 < c do i=1,ndih_constr < do i=idihconstr_start,idihconstr_end < itori=idih_constr(i) < phii=phi(itori) < difi=pinorm(phii-phi0(i)) < if (difi.gt.drange(i)) then < difi=difi-drange(i) < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 < else if (difi.lt.-drange(i)) then < difi=difi+drange(i) < edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 < gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 < else < difi=0.0 < endif < enddo < endif < return < end < < C The rigorous attempt to derive energy function < subroutine ebend_kcc(etheta,ethetacnstr) < < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.VAR' < include 'COMMON.GEO' < include 'COMMON.LOCAL' < include 'COMMON.TORSION' < include 'COMMON.INTERACT' < include 'COMMON.DERIV' < include 'COMMON.CHAIN' < include 'COMMON.NAMES' < include 'COMMON.IOUNITS' < include 'COMMON.FFIELD' < include 'COMMON.TORCNSTR' < include 'COMMON.CONTROL' < logical lprn < double precision thybt1(maxtermkcc) < C Set lprn=.true. for debugging < lprn=.false. < c lprn=.true. < C print *,"wchodze kcc" < if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode < if (tor_mode.ne.2) etheta=0.0D0 < do i=ithet_start,ithet_end < c print *,i,itype(i-1),itype(i),itype(i-2) < if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 < & .or.itype(i).eq.ntyp1) cycle < iti=itortyp_kcc(itype(i-1)) < sinthet=dsin(theta(i)/2.0d0) < costhet=dcos(theta(i)/2.0d0) < do j=1,nbend_kcc_Tb(iti) < thybt1(j)=v1bend_chyb(j,iti) < enddo < sumth1thyb=tschebyshev < & (1,nbend_kcc_Tb(iti),thybt1(1),costhet) < if (lprn) write (iout,*) i-1,itype(i-1),iti,theta(i)*rad2deg, < & sumth1thyb < ihelp=nbend_kcc_Tb(iti)-1 < gradthybt1=gradtschebyshev < & (0,ihelp,thybt1(1),costhet) < etheta=etheta+sumth1thyb < C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0) < gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang* < & gradthybt1*sinthet*(-0.5d0) < enddo < if (tor_mode.ne.2) then < ethetacnstr=0.0d0 < C print *,ithetaconstr_start,ithetaconstr_end,"TU" < do i=ithetaconstr_start,ithetaconstr_end < itheta=itheta_constr(i) < thetiii=theta(itheta) < difi=pinorm(thetiii-theta_constr0(i)) < if (difi.gt.theta_drange(i)) then < difi=difi-theta_drange(i) < ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 < gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) < & +for_thet_constr(i)*difi**3 < else if (difi.lt.-drange(i)) then < difi=difi+drange(i) < ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 < gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) < & +for_thet_constr(i)*difi**3 < else < difi=0.0 < endif < if (energy_dec) then < write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", < & i,itheta,rad2deg*thetiii, < & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), < & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, < & gloc(itheta+nphi-2,icg) < endif < enddo < endif < return < end 7845d6103 < include 'COMMON.SHIELD' 8264d6521 < include 'COMMON.SHIELD' 8567d6823 < CC & *fac_shield(i)**2*fac_shield(j)**2 8687,8688d6942 < include 'COMMON.SHIELD' < include 'COMMON.CONTROL' 8692d6945 < C print *,"wchodze",fac_shield(i),shield_mode 8701,8702d6953 < C* < C & fac_shield(i)**2*fac_shield(j)**2 8766d7016 < C print *,ekont,ees,i,k 8768,8848d7017 < C now gradient over shielding < C return < if (shield_mode.gt.0) then < j=ees0plist(jj,i) < l=ees0plist(kk,k) < C print *,i,j,fac_shield(i),fac_shield(j), < C &fac_shield(k),fac_shield(l) < if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. < & (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then < do ilist=1,ishield_list(i) < iresshield=shield_list(ilist,i) < do m=1,3 < rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i) < C & *2.0 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ < & rlocshield < & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i) < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) < &+rlocshield < enddo < enddo < do ilist=1,ishield_list(j) < iresshield=shield_list(ilist,j) < do m=1,3 < rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j) < C & *2.0 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ < & rlocshield < & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j) < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) < & +rlocshield < enddo < enddo < < do ilist=1,ishield_list(k) < iresshield=shield_list(ilist,k) < do m=1,3 < rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k) < C & *2.0 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ < & rlocshield < & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k) < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) < & +rlocshield < enddo < enddo < do ilist=1,ishield_list(l) < iresshield=shield_list(ilist,l) < do m=1,3 < rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l) < C & *2.0 < gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+ < & rlocshield < & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l) < gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1) < & +rlocshield < enddo < enddo < C print *,gshieldx(m,iresshield) < do m=1,3 < gshieldc_ec(m,i)=gshieldc_ec(m,i)+ < & grad_shield(m,i)*ehbcorr/fac_shield(i) < gshieldc_ec(m,j)=gshieldc_ec(m,j)+ < & grad_shield(m,j)*ehbcorr/fac_shield(j) < gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+ < & grad_shield(m,i)*ehbcorr/fac_shield(i) < gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+ < & grad_shield(m,j)*ehbcorr/fac_shield(j) < < gshieldc_ec(m,k)=gshieldc_ec(m,k)+ < & grad_shield(m,k)*ehbcorr/fac_shield(k) < gshieldc_ec(m,l)=gshieldc_ec(m,l)+ < & grad_shield(m,l)*ehbcorr/fac_shield(l) < gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+ < & grad_shield(m,k)*ehbcorr/fac_shield(k) < gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+ < & grad_shield(m,l)*ehbcorr/fac_shield(l) < < enddo < endif < endif 8867c7036 < iti1 = itortyp(itype(i+1)) --- > iti1 = itype2loc(itype(i+1)) 9531c7700 < call transpose2(EE(1,1,k),auxmat(1,1)) --- > call transpose2(EE(1,1,itk),auxmat(1,1)) 9612c7781 < call transpose2(EE(1,1,l),auxmat(1,1)) --- > call transpose2(EE(1,1,itl),auxmat(1,1)) 9685c7854 < call transpose2(EE(1,1,j),auxmat(1,1)) --- > call transpose2(EE(1,1,itj),auxmat(1,1)) 9773c7942 < gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2) --- > gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2) 9775c7944 < gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2) --- > gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2) 10080c8249 < C o| o | | o |o C --- > C o| o | | o |o C 10084,10085c8253,8254 < C i i C < C C --- > C i i C > C C 10256c8425 < C C --- > C C 10259c8428 < C o o C --- > C o o C 10272c8441 < iti=itortyp(itype(i)) --- > iti=itype2loc(itype(i)) 10292c8461 < call transpose2(EE(1,1,k),auxmat(1,1)) --- > call transpose2(EE(1,1,itk),auxmat(1,1)) 10373c8542 < C C --- > C C 10381c8550 < C \ / \ \ / \ C --- > C \ / \ \ / \ C 10384c8553 < C C --- > C C 11082,11656d9250 < return < end < CCC---------------------------------------------- < subroutine Eliptransfer(eliptran) < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.GEO' < include 'COMMON.VAR' < include 'COMMON.LOCAL' < include 'COMMON.CHAIN' < include 'COMMON.DERIV' < include 'COMMON.NAMES' < include 'COMMON.INTERACT' < include 'COMMON.IOUNITS' < include 'COMMON.CALC' < include 'COMMON.CONTROL' < include 'COMMON.SPLITELE' < include 'COMMON.SBRIDGE' < C this is done by Adasko < C print *,"wchodze" < C structure of box: < C water < C--bordliptop-- buffore starts < C--bufliptop--- here true lipid starts < C lipid < C--buflipbot--- lipid ends buffore starts < C--bordlipbot--buffore ends < eliptran=0.0 < do i=ilip_start,ilip_end < C do i=1,1 < if (itype(i).eq.ntyp1) cycle < < positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) < if (positi.le.0) positi=positi+boxzsize < C print *,i < C first for peptide groups < c for each residue check if it is in lipid or lipid water border area < if ((positi.gt.bordlipbot) < &.and.(positi.lt.bordliptop)) then < C the energy transfer exist < if (positi.lt.buflipbot) then < C what fraction I am in < fracinbuf=1.0d0- < & ((positi-bordlipbot)/lipbufthick) < C lipbufthick is thickenes of lipid buffore < sslip=sscalelip(fracinbuf) < ssgradlip=-sscagradlip(fracinbuf)/lipbufthick < eliptran=eliptran+sslip*pepliptran < gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 < gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 < C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran < < C print *,"doing sccale for lower part" < C print *,i,sslip,fracinbuf,ssgradlip < elseif (positi.gt.bufliptop) then < fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) < sslip=sscalelip(fracinbuf) < ssgradlip=sscagradlip(fracinbuf)/lipbufthick < eliptran=eliptran+sslip*pepliptran < gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0 < gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 < C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran < C print *, "doing sscalefor top part" < C print *,i,sslip,fracinbuf,ssgradlip < else < eliptran=eliptran+pepliptran < C print *,"I am in true lipid" < endif < C else < C eliptran=elpitran+0.0 ! I am in water < endif < enddo < C print *, "nic nie bylo w lipidzie?" < C now multiply all by the peptide group transfer factor < C eliptran=eliptran*pepliptran < C now the same for side chains < CV do i=1,1 < do i=ilip_start,ilip_end < if (itype(i).eq.ntyp1) cycle < positi=(mod(c(3,i+nres),boxzsize)) < if (positi.le.0) positi=positi+boxzsize < C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop < c for each residue check if it is in lipid or lipid water border area < C respos=mod(c(3,i+nres),boxzsize) < C print *,positi,bordlipbot,buflipbot < if ((positi.gt.bordlipbot) < & .and.(positi.lt.bordliptop)) then < C the energy transfer exist < if (positi.lt.buflipbot) then < fracinbuf=1.0d0- < & ((positi-bordlipbot)/lipbufthick) < C lipbufthick is thickenes of lipid buffore < sslip=sscalelip(fracinbuf) < ssgradlip=-sscagradlip(fracinbuf)/lipbufthick < eliptran=eliptran+sslip*liptranene(itype(i)) < gliptranx(3,i)=gliptranx(3,i) < &+ssgradlip*liptranene(itype(i)) < gliptranc(3,i-1)= gliptranc(3,i-1) < &+ssgradlip*liptranene(itype(i)) < C print *,"doing sccale for lower part" < elseif (positi.gt.bufliptop) then < fracinbuf=1.0d0- < &((bordliptop-positi)/lipbufthick) < sslip=sscalelip(fracinbuf) < ssgradlip=sscagradlip(fracinbuf)/lipbufthick < eliptran=eliptran+sslip*liptranene(itype(i)) < gliptranx(3,i)=gliptranx(3,i) < &+ssgradlip*liptranene(itype(i)) < gliptranc(3,i-1)= gliptranc(3,i-1) < &+ssgradlip*liptranene(itype(i)) < C print *, "doing sscalefor top part",sslip,fracinbuf < else < eliptran=eliptran+liptranene(itype(i)) < C print *,"I am in true lipid" < endif < endif ! if in lipid or buffor < C else < C eliptran=elpitran+0.0 ! I am in water < enddo < return < end < C--------------------------------------------------------- < C AFM soubroutine for constant force < subroutine AFMforce(Eafmforce) < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.GEO' < include 'COMMON.VAR' < include 'COMMON.LOCAL' < include 'COMMON.CHAIN' < include 'COMMON.DERIV' < include 'COMMON.NAMES' < include 'COMMON.INTERACT' < include 'COMMON.IOUNITS' < include 'COMMON.CALC' < include 'COMMON.CONTROL' < include 'COMMON.SPLITELE' < include 'COMMON.SBRIDGE' < real*8 diffafm(3) < dist=0.0d0 < Eafmforce=0.0d0 < do i=1,3 < diffafm(i)=c(i,afmend)-c(i,afmbeg) < dist=dist+diffafm(i)**2 < enddo < dist=dsqrt(dist) < Eafmforce=-forceAFMconst*(dist-distafminit) < do i=1,3 < gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist < gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist < enddo < C print *,'AFM',Eafmforce < return < end < C--------------------------------------------------------- < C AFM subroutine with pseudoconstant velocity < subroutine AFMvel(Eafmforce) < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.GEO' < include 'COMMON.VAR' < include 'COMMON.LOCAL' < include 'COMMON.CHAIN' < include 'COMMON.DERIV' < include 'COMMON.NAMES' < include 'COMMON.INTERACT' < include 'COMMON.IOUNITS' < include 'COMMON.CALC' < include 'COMMON.CONTROL' < include 'COMMON.SPLITELE' < include 'COMMON.SBRIDGE' < real*8 diffafm(3) < C Only for check grad COMMENT if not used for checkgrad < C totT=3.0d0 < C-------------------------------------------------------- < C print *,"wchodze" < dist=0.0d0 < Eafmforce=0.0d0 < do i=1,3 < diffafm(i)=c(i,afmend)-c(i,afmbeg) < dist=dist+diffafm(i)**2 < enddo < dist=dsqrt(dist) < Eafmforce=0.5d0*forceAFMconst < & *(distafminit+totTafm*velAFMconst-dist)**2 < C Eafmforce=-forceAFMconst*(dist-distafminit) < do i=1,3 < gradafm(i,afmend-1)=-forceAFMconst* < &(distafminit+totTafm*velAFMconst-dist) < &*diffafm(i)/dist < gradafm(i,afmbeg-1)=forceAFMconst* < &(distafminit+totTafm*velAFMconst-dist) < &*diffafm(i)/dist < enddo < C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist < return < end < C----------------------------------------------------------- < C first for shielding is setting of function of side-chains < subroutine set_shield_fac < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.CHAIN' < include 'COMMON.DERIV' < include 'COMMON.IOUNITS' < include 'COMMON.SHIELD' < include 'COMMON.INTERACT' < C this is the squar root 77 devided by 81 the epislion in lipid (in protein) < double precision div77_81/0.974996043d0/, < &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) < < C the vector between center of side_chain and peptide group < double precision pep_side(3),long,side_calf(3), < &pept_group(3),costhet_grad(3),cosphi_grad_long(3), < &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) < C the line belowe needs to be changed for FGPROC>1 < do i=1,nres-1 < if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle < ishield_list(i)=0 < Cif there two consequtive dummy atoms there is no peptide group between them < C the line below has to be changed for FGPROC>1 < VolumeTotal=0.0 < do k=1,nres < if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle < dist_pep_side=0.0 < dist_side_calf=0.0 < do j=1,3 < C first lets set vector conecting the ithe side-chain with kth side-chain < pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 < C pep_side(j)=2.0d0 < C and vector conecting the side-chain with its proper calfa < side_calf(j)=c(j,k+nres)-c(j,k) < C side_calf(j)=2.0d0 < pept_group(j)=c(j,i)-c(j,i+1) < C lets have their lenght < dist_pep_side=pep_side(j)**2+dist_pep_side < dist_side_calf=dist_side_calf+side_calf(j)**2 < dist_pept_group=dist_pept_group+pept_group(j)**2 < enddo < dist_pep_side=dsqrt(dist_pep_side) < dist_pept_group=dsqrt(dist_pept_group) < dist_side_calf=dsqrt(dist_side_calf) < do j=1,3 < pep_side_norm(j)=pep_side(j)/dist_pep_side < side_calf_norm(j)=dist_side_calf < enddo < C now sscale fraction < sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield < C print *,buff_shield,"buff" < C now sscale < if (sh_frac_dist.le.0.0) cycle < C If we reach here it means that this side chain reaches the shielding sphere < C Lets add him to the list for gradient < ishield_list(i)=ishield_list(i)+1 < C ishield_list is a list of non 0 side-chain that contribute to factor gradient < C this list is essential otherwise problem would be O3 < shield_list(ishield_list(i),i)=k < C Lets have the sscale value < if (sh_frac_dist.gt.1.0) then < scale_fac_dist=1.0d0 < do j=1,3 < sh_frac_dist_grad(j)=0.0d0 < enddo < else < scale_fac_dist=-sh_frac_dist*sh_frac_dist < & *(2.0*sh_frac_dist-3.0d0) < fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2) < & /dist_pep_side/buff_shield*0.5 < C remember for the final gradient multiply sh_frac_dist_grad(j) < C for side_chain by factor -2 ! < do j=1,3 < sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) < C print *,"jestem",scale_fac_dist,fac_help_scale, < C & sh_frac_dist_grad(j) < enddo < endif < C if ((i.eq.3).and.(k.eq.2)) then < C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist < C & ,"TU" < C endif < < C this is what is now we have the distance scaling now volume... < short=short_r_sidechain(itype(k)) < long=long_r_sidechain(itype(k)) < costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2) < C now costhet_grad < C costhet=0.0d0 < costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4 < C costhet_fac=0.0d0 < do j=1,3 < costhet_grad(j)=costhet_fac*pep_side(j) < enddo < C remember for the final gradient multiply costhet_grad(j) < C for side_chain by factor -2 ! < C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 < C pep_side0pept_group is vector multiplication < pep_side0pept_group=0.0 < do j=1,3 < pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) < enddo < cosalfa=(pep_side0pept_group/ < & (dist_pep_side*dist_side_calf)) < fac_alfa_sin=1.0-cosalfa**2 < fac_alfa_sin=dsqrt(fac_alfa_sin) < rkprim=fac_alfa_sin*(long-short)+short < C now costhet_grad < cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2) < cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4 < < do j=1,3 < cosphi_grad_long(j)=cosphi_fac*pep_side(j) < &+cosphi**3*0.5/dist_pep_side**2*(-rkprim) < &*(long-short)/fac_alfa_sin*cosalfa/ < &((dist_pep_side*dist_side_calf))* < &((side_calf(j))-cosalfa* < &((pep_side(j)/dist_pep_side)*dist_side_calf)) < < cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim) < &*(long-short)/fac_alfa_sin*cosalfa < &/((dist_pep_side*dist_side_calf))* < &(pep_side(j)- < &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) < enddo < < VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi) < & /VSolvSphere_div < & *wshield < C now the gradient... < C grad_shield is gradient of Calfa for peptide groups < C write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist, < C & costhet,cosphi < C write(iout,*) "cosphi_compon",i,k,pep_side0pept_group, < C & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k) < do j=1,3 < grad_shield(j,i)=grad_shield(j,i) < C gradient po skalowaniu < & +(sh_frac_dist_grad(j) < C gradient po costhet < &-scale_fac_dist*costhet_grad(j)/(1.0-costhet) < &-scale_fac_dist*(cosphi_grad_long(j)) < &/(1.0-cosphi) )*div77_81 < &*VofOverlap < C grad_shield_side is Cbeta sidechain gradient < grad_shield_side(j,ishield_list(i),i)= < & (sh_frac_dist_grad(j)*-2.0d0 < & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet) < & +scale_fac_dist*(cosphi_grad_long(j)) < & *2.0d0/(1.0-cosphi)) < & *div77_81*VofOverlap < < grad_shield_loc(j,ishield_list(i),i)= < & scale_fac_dist*cosphi_grad_loc(j) < & *2.0d0/(1.0-cosphi) < & *div77_81*VofOverlap < enddo < VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist < enddo < fac_shield(i)=VolumeTotal*div77_81+div4_81 < C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) < enddo < return < end < C-------------------------------------------------------------------------- < double precision function tschebyshev(m,n,x,y) < implicit none < include "DIMENSIONS" < integer i,m,n < double precision x(n),y,yy(0:maxvar),aux < c Tschebyshev polynomial. Note that the first term is omitted < c m=0: the constant term is included < c m=1: the constant term is not included < yy(0)=1.0d0 < yy(1)=y < do i=2,n < yy(i)=2*yy(1)*yy(i-1)-yy(i-2) < enddo < aux=0.0d0 < do i=m,n < aux=aux+x(i)*yy(i) < enddo < tschebyshev=aux < return < end < C-------------------------------------------------------------------------- < double precision function gradtschebyshev(m,n,x,y) < implicit none < include "DIMENSIONS" < integer i,m,n < double precision x(n+1),y,yy(0:maxvar),aux < c Tschebyshev polynomial. Note that the first term is omitted < c m=0: the constant term is included < c m=1: the constant term is not included < yy(0)=1.0d0 < yy(1)=2.0d0*y < do i=2,n < yy(i)=2*y*yy(i-1)-yy(i-2) < enddo < aux=0.0d0 < do i=m,n < aux=aux+x(i+1)*yy(i)*(i+1) < C print *, x(i+1),yy(i),i < enddo < gradtschebyshev=aux < return < end < C------------------------------------------------------------------------ < C first for shielding is setting of function of side-chains < subroutine set_shield_fac2 < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < include 'COMMON.CHAIN' < include 'COMMON.DERIV' < include 'COMMON.IOUNITS' < include 'COMMON.SHIELD' < include 'COMMON.INTERACT' < C this is the squar root 77 devided by 81 the epislion in lipid (in protein) < double precision div77_81/0.974996043d0/, < &div4_81/0.2222222222d0/,sh_frac_dist_grad(3) < < C the vector between center of side_chain and peptide group < double precision pep_side(3),long,side_calf(3), < &pept_group(3),costhet_grad(3),cosphi_grad_long(3), < &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) < C the line belowe needs to be changed for FGPROC>1 < do i=1,nres-1 < if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle < ishield_list(i)=0 < Cif there two consequtive dummy atoms there is no peptide group between them < C the line below has to be changed for FGPROC>1 < VolumeTotal=0.0 < do k=1,nres < if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle < dist_pep_side=0.0 < dist_side_calf=0.0 < do j=1,3 < C first lets set vector conecting the ithe side-chain with kth side-chain < pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0 < C pep_side(j)=2.0d0 < C and vector conecting the side-chain with its proper calfa < side_calf(j)=c(j,k+nres)-c(j,k) < C side_calf(j)=2.0d0 < pept_group(j)=c(j,i)-c(j,i+1) < C lets have their lenght < dist_pep_side=pep_side(j)**2+dist_pep_side < dist_side_calf=dist_side_calf+side_calf(j)**2 < dist_pept_group=dist_pept_group+pept_group(j)**2 < enddo < dist_pep_side=dsqrt(dist_pep_side) < dist_pept_group=dsqrt(dist_pept_group) < dist_side_calf=dsqrt(dist_side_calf) < do j=1,3 < pep_side_norm(j)=pep_side(j)/dist_pep_side < side_calf_norm(j)=dist_side_calf < enddo < C now sscale fraction < sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield < C print *,buff_shield,"buff" < C now sscale < if (sh_frac_dist.le.0.0) cycle < C If we reach here it means that this side chain reaches the shielding sphere < C Lets add him to the list for gradient < ishield_list(i)=ishield_list(i)+1 < C ishield_list is a list of non 0 side-chain that contribute to factor gradient < C this list is essential otherwise problem would be O3 < shield_list(ishield_list(i),i)=k < C Lets have the sscale value < if (sh_frac_dist.gt.1.0) then < scale_fac_dist=1.0d0 < do j=1,3 < sh_frac_dist_grad(j)=0.0d0 < enddo < else < scale_fac_dist=-sh_frac_dist*sh_frac_dist < & *(2.0d0*sh_frac_dist-3.0d0) < fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2) < & /dist_pep_side/buff_shield*0.5d0 < C remember for the final gradient multiply sh_frac_dist_grad(j) < C for side_chain by factor -2 ! < do j=1,3 < sh_frac_dist_grad(j)=fac_help_scale*pep_side(j) < C sh_frac_dist_grad(j)=0.0d0 < C scale_fac_dist=1.0d0 < C print *,"jestem",scale_fac_dist,fac_help_scale, < C & sh_frac_dist_grad(j) < enddo < endif < C this is what is now we have the distance scaling now volume... < short=short_r_sidechain(itype(k)) < long=long_r_sidechain(itype(k)) < costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) < sinthet=short/dist_pep_side*costhet < C now costhet_grad < C costhet=0.6d0 < C sinthet=0.8 < costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4 < C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet < C & -short/dist_pep_side**2/costhet) < C costhet_fac=0.0d0 < do j=1,3 < costhet_grad(j)=costhet_fac*pep_side(j) < enddo < C remember for the final gradient multiply costhet_grad(j) < C for side_chain by factor -2 ! < C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1 < C pep_side0pept_group is vector multiplication < pep_side0pept_group=0.0d0 < do j=1,3 < pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) < enddo < cosalfa=(pep_side0pept_group/ < & (dist_pep_side*dist_side_calf)) < fac_alfa_sin=1.0d0-cosalfa**2 < fac_alfa_sin=dsqrt(fac_alfa_sin) < rkprim=fac_alfa_sin*(long-short)+short < C rkprim=short < < C now costhet_grad < cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2) < C cosphi=0.6 < cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4 < sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/ < & dist_pep_side**2) < C sinphi=0.8 < do j=1,3 < cosphi_grad_long(j)=cosphi_fac*pep_side(j) < &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) < &*(long-short)/fac_alfa_sin*cosalfa/ < &((dist_pep_side*dist_side_calf))* < &((side_calf(j))-cosalfa* < &((pep_side(j)/dist_pep_side)*dist_side_calf)) < C cosphi_grad_long(j)=0.0d0 < cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) < &*(long-short)/fac_alfa_sin*cosalfa < &/((dist_pep_side*dist_side_calf))* < &(pep_side(j)- < &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side) < C cosphi_grad_loc(j)=0.0d0 < enddo < C print *,sinphi,sinthet < VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) < & /VSolvSphere_div < C & *wshield < C now the gradient... < do j=1,3 < grad_shield(j,i)=grad_shield(j,i) < C gradient po skalowaniu < & +(sh_frac_dist_grad(j)*VofOverlap < C gradient po costhet < & +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0* < &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( < & sinphi/sinthet*costhet*costhet_grad(j) < & +sinthet/sinphi*cosphi*cosphi_grad_long(j))) < & )*wshield < C grad_shield_side is Cbeta sidechain gradient < grad_shield_side(j,ishield_list(i),i)= < & (sh_frac_dist_grad(j)*-2.0d0 < & *VofOverlap < & -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0* < &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( < & sinphi/sinthet*costhet*costhet_grad(j) < & +sinthet/sinphi*cosphi*cosphi_grad_long(j))) < & )*wshield < < grad_shield_loc(j,ishield_list(i),i)= < & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0* < &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*( < & sinthet/sinphi*cosphi*cosphi_grad_loc(j) < & )) < & *wshield < enddo < VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist < enddo < fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) < C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) < enddo