X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=c90fc3ac271ba7bdb909db59a48007dac278b292;hb=6c0d9d49feaa4dd54945a4e74fa72a18efdf3a0e;hp=49676d8b767988768b851732bbc2c35e1c47c13d;hpb=e910def1e1c1c0f0b35e8ce1cbf2a9ee59628771;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 49676d8..c90fc3a 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -141,9 +141,11 @@ 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 - write (iout,*) "shield_mode",shield_mode - if (shield_mode.gt.0) then +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 c print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then @@ -201,7 +203,14 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then + 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 else ebe=0 ethetacnstr=0 @@ -217,8 +226,16 @@ C C Calculate the virtual-bond torsional energy. C cd print *,'nterm=',nterm +C print *,"tor",tor_mode if (wtor.gt.0) then + if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then call etor(etors,edihcnstr) + 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 else etors=0 edihcnstr=0 @@ -227,7 +244,7 @@ c print *,"Processor",myrank," computed Utor" C C 6/23/01 Calculate double-torsional energy C - if (wtor_d.gt.0) then + if ((wtor_d.gt.0).and.((tor_mode.eq.0).or.(tor_mode.eq.2))) then call etor_d(etors_d) else etors_d=0 @@ -546,6 +563,12 @@ c enddo & wstrain*ghpbc(j,i) & +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) + enddo enddo @@ -564,6 +587,11 @@ c enddo & wstrain*ghpbc(j,i) & +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) + enddo enddo @@ -682,6 +710,13 @@ c enddo do i=-1,nct do j=1,3 #ifdef SPLITELE +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" gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ @@ -702,11 +737,27 @@ c enddo & +wscloc*gscloc(j,i) & +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) + + + + + + #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ & 0.5d0*(wscp*gvdwc_scpp(j,i)+ - & welec*gelc_long(j,i) + & welec*gelc_long(j,i)+ & wel_loc*gel_loc_long(j,i)+ & wcorr*gcorr_long(j,i)+ & wcorr5*gradcorr5_long(j,i)+ @@ -723,6 +774,20 @@ c enddo & +wscloc*gscloc(j,i) & +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) + + + + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ @@ -731,6 +796,14 @@ c enddo & wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) & +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) + + + enddo enddo #ifdef DEBUG @@ -933,6 +1006,7 @@ c------------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' + include 'COMMON.CONTROL' double precision kfac /2.4d0/ double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/ c facT=temp0/t_bath @@ -968,6 +1042,11 @@ c facT=2*temp0/(t_bath+temp0) #endif stop 555 endif + if (shield_mode.gt.0) then + wscp=weights(2)*fact + wsc=weights(1)*fact + wvdwpp=weights(16)*fact + endif welec=weights(3)*fact wcorr=weights(4)*fact3 wcorr5=weights(5)*fact4 @@ -1028,7 +1107,7 @@ C------------------------------------------------------------------------ & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, - & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, + & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr, & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, & etot 10 format (/'Virtual-chain energies:'// @@ -2704,28 +2783,28 @@ c write(iout,*) 'nphi=',nphi,nres #endif #ifdef NEWCORR if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp+1 + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp+1 + iti1=nloctyp endif c write(iout,*),i - b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0) + b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0) & +bnew1(2,1,iti)*dsin(theta(i-1)) - & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0) + & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0) gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 & +bnew1(2,1,iti)*dcos(theta(i-1)) & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0) + b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0) & +bnew2(2,1,iti)*dsin(theta(i-1)) - & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0) + & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0) c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 @@ -2763,6 +2842,17 @@ c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) enddo #else + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itype2loc(itype(i-2)) + else + iti=nloctyp + endif +c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itype2loc(itype(i-1)) + else + iti1=nloctyp + endif b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) b2(1,i-2)=b(2,iti) @@ -2850,15 +2940,15 @@ c write (iout,*) 'theta=', theta(i-1) endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp + iti1=nloctyp endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) @@ -2871,8 +2961,8 @@ c if (i .gt. iatel_s+2) then call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" #endif -c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti), -c & EE(1,2,iti),EE(2,2,iti) +c write(iout,*) "co jest kurwa", iti, EE(1,1,i),EE(2,1,i), +c & EE(1,2,iti),EE(2,2,i) call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) c write(iout,*) "Macierz EUG", @@ -2907,17 +2997,24 @@ c & eug(2,2,i-2) c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then if (itype(i-1).le.ntyp) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp + iti1=nloctyp endif else - iti1=ntortyp + iti1=nloctyp endif do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -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,*) 'mu1',mu1(:,i-2) cd write (iout,*) 'mu2',mu2(:,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) @@ -3204,11 +3301,11 @@ c endif #endif #endif cd do i=1,nres -cd iti = itortyp(itype(i)) +cd iti = itype2loc(itype(i)) cd write (iout,*) i cd do j=1,2 cd write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)') -cd & (EE(j,k,iti),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2) +cd & (EE(j,k,i),k=1,2),(Ug(j,k,i),k=1,2),(EUg(j,k,i),k=1,2) cd enddo cd enddo return @@ -3326,21 +3423,23 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end - if (i.le.1) cycle +c if (i.le.1) cycle C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+4).gt.nres) - & .or.((i-1).le.0) +C Adam: Unnecessary: handled by iturn3_end and iturn3_start +c & .or.((i+4).gt.nres) +c & .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 +C Adam: Instructions below will switch off existing interactions +c if(i.gt.1)then +c if(itype(i-1).eq.ntyp1)cycle +c end if +c if(i.LT.nres-3)then +c if (itype(i+4).eq.ntyp1) cycle +c end if dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3365,14 +3464,14 @@ C end of changes by Ana 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+5).gt.nres) - & .or.((i-1).le.0) +c & .or.((i+5).gt.nres) +c & .or.((i-1).le.0) C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +c & .or. itype(i+5).eq.ntyp1 +c & .or. itype(i).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3429,15 +3528,17 @@ C do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c +CTU KURWA do i=iatel_s,iatel_e - if (i.le.1) cycle +C do i=75,75 +c 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 & .or.((i+2).gt.nres) +c & .or.((i-1).le.0) C end of changes by Ana - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +c & .or. itype(i+2).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3486,16 +3587,18 @@ c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) +C I TU KURWA do j=ielstart(i),ielend(i) +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 & .or.((j+2).gt.nres) +c & .or.((j-1).le.0) C end of changes by Ana - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +c & .or.itype(j+2).eq.ntyp1 +c & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j @@ -3669,12 +3772,20 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac C MARYSIA - eesij=(el1+el2) +C eesij=(el1+el2) C 12/26/95 - for the evaluation of multi-body H-bonding interactions ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg) if (shield_mode.gt.0) then - ees=ees+eesij*fac_shield(i)*fac_shield(j) +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) ees=ees+eesij endif evdw1=evdw1+evdwij*sss @@ -3687,7 +3798,8 @@ cd & xmedi,ymedi,zmedi,xj,yj,zj write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &'evdw1',i,j,evdwij &,iteli,itelj,aaa,evdw1 - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij + write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, + &fac_shield(i),fac_shield(j) endif C @@ -3700,22 +3812,101 @@ C erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij + * * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=facel*xj ggg(2)=facel*yj ggg(3)=facel*zj + 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 c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf c gelc(k,j)=gelc(k,j)+ghalf c enddo c 9/28/08 AL Gradient compotents will be summed only at the end +C print *,"before", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc_long(k,j)=gelc_long(k,j)+ggg(k) +C & +grad_shield(k,j)*eesij/fac_shield(j) gelc_long(k,i)=gelc_long(k,i)-ggg(k) +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) enddo +C print *,"bafter", gelc_long(1,i), gelc_long(1,j) + * * Loop over residues i+1 thru j-1. * @@ -3764,8 +3955,11 @@ C MARYSIA * Radial derivatives. First process both termini of the fragment (i,j) * ggg(1)=fac*xj +C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j) ggg(2)=fac*yj +C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j) ggg(3)=fac*zj +C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j) c do k=1,3 c ghalf=0.5D0*ggg(k) c gelc(k,i)=gelc(k,i)+ghalf @@ -3808,7 +4002,8 @@ c 9/28/08 AL Gradient compotents will be summed only at the end cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 - ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) + ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* + & fac_shield(i)**2*fac_shield(j)**2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -3824,16 +4019,21 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo +C print *,"before22", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc(k,i)=gelc(k,i) - & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) + & *fac_shield(i)**2*fac_shield(j)**2 gelc(k,j)=gelc(k,j) - & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) + & *fac_shield(i)**2*fac_shield(j)**2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo +C print *,"before33", gelc_long(1,i), gelc_long(1,j) + C MARYSIA c endif !sscale IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 @@ -4034,15 +4234,71 @@ cgrad endif C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) + 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=',muij(1),muij(2),muij(3),muij(4) +C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4) C Calculate patrial derivative for theta angle #ifdef NEWCORR - geel_loc_ij=a22*gmuij1(1) + geel_loc_ij=(a22*gmuij1(1) & +a23*gmuij1(2) & +a32*gmuij1(3) - & +a33*gmuij1(4) + & +a33*gmuij1(4)) + & *fac_shield(i)*fac_shield(j) c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4058,6 +4314,8 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc + & *fac_shield(i)*fac_shield(j) + c Derivative over j residue geel_loc_ji=a22*gmuji1(1) & +a23*gmuji1(2) @@ -4069,6 +4327,8 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc + & *fac_shield(i)*fac_shield(j) + geel_loc_ji= & +a22*gmuji2(1) & +a23*gmuji2(2) @@ -4079,6 +4339,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc + & *fac_shield(i)*fac_shield(j) #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4092,15 +4353,19 @@ c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) C Partial derivatives in virtual-bond dihedral angles gamma if (i.gt.1) & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ - & 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) + & (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) + gel_loc_loc(j-1)=gel_loc_loc(j-1)+ - & 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) + & (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) C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 - ggg(l)=agg(l,1)*muij(1)+ - & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4) + 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) gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) @@ -4116,12 +4381,20 @@ C Remaining derivatives of eello do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) + & *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) + enddo ENDIF C Change 12/26/95 to calculate four-body contributions to H-bonding energy @@ -4200,8 +4473,18 @@ c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) ees0mij=0 endif c ees0mij=0.0D0 + 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 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) + & *fac_shield(i)*fac_shield(j) ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) + & *fac_shield(i)*fac_shield(j) C Diagnostics. Comment out or remove after debugging! c ees0p(num_conti,i)=0.5D0*fac3*ees0pij c ees0m(num_conti,i)=0.5D0*fac3*ees0mij @@ -4269,17 +4552,29 @@ cgrad ghalfm=0.5D0*gggm(k) gacontp_hb1(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *fac_shield(i)*fac_shield(j) + gacontp_hb2(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *fac_shield(i)*fac_shield(j) + gacontp_hb3(k,num_conti,i)=gggp(k) + & *fac_shield(i)*fac_shield(j) + gacontm_hb1(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) + & *fac_shield(i)*fac_shield(j) + gacontm_hb2(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) + & *fac_shield(i)*fac_shield(j) + gacontm_hb3(k,num_conti,i)=gggm(k) + & *fac_shield(i)*fac_shield(j) + enddo C Diagnostics. Comment out or remove after debugging! cdiag do k=1,3 @@ -4331,6 +4626,7 @@ C Third- and fourth-order contributions from turns include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.SHIELD' dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), @@ -4371,15 +4667,71 @@ c auxalary matrix for i+2 and constant i+1 call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1)) call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1)) + 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 eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) C Derivatives in theta gloc(nphi+i,icg)=gloc(nphi+i,icg) & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 + & *fac_shield(i)*fac_shield(j) gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 + & *fac_shield(i)*fac_shield(j) - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) + +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 + + 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)') cd write (2,*) 'i,',i,' j',j,'eello_turn3', cd & 0.5d0*(pizda(1,1)+pizda(2,2)), cd & ' eello_turn3_num',4*eello_turn3_num @@ -4388,12 +4740,14 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),auxmat3(1,1)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -4407,6 +4761,8 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) + a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) a_temp(2,1)=aggi1(l,3)!+agg(l,3) @@ -4414,6 +4770,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -4421,6 +4778,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -4428,6 +4786,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) + & *fac_shield(i)*fac_shield(j) enddo return end @@ -4448,6 +4807,7 @@ C Third- and fourth-order contributions from turns include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.SHIELD' dimension ggg(3) double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), @@ -4480,9 +4840,9 @@ c write(iout,*)"WCHODZE W PROGRAM" a_temp(1,2)=a23 a_temp(2,1)=a32 a_temp(2,2)=a33 - iti1=itortyp(itype(i+1)) - iti2=itortyp(itype(i+2)) - iti3=itortyp(itype(i+3)) + iti1=itype2loc(itype(i+1)) + iti2=itype2loc(itype(i+2)) + iti3=itype2loc(itype(i+3)) c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) @@ -4548,20 +4908,82 @@ c i+3 gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2)) gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2)) gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2)) - + 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 + + + + + + cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), cd & ' eello_turn4_num',8*eello_turn4_num #ifdef NEWCORR gloc(nphi+i,icg)=gloc(nphi+i,icg) & -(gs13+gsE13+gsEE1)*wturn4 + & *fac_shield(i)*fac_shield(j) gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg) & -(gs23+gs21+gsEE2)*wturn4 + & *fac_shield(i)*fac_shield(j) + gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg) & -(gs32+gsE31+gsEE3)*wturn4 + & *fac_shield(i)*fac_shield(j) + c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- c & gs2 #endif @@ -4577,6 +4999,7 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) + & *fac_shield(i)*fac_shield(j) C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -4585,6 +5008,7 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) + & *fac_shield(i)*fac_shield(j) C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -4596,6 +5020,7 @@ C Derivatives in gamma(i+2) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -4615,6 +5040,7 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) enddo endif C Remaining derivatives of this turn contribution @@ -4633,6 +5059,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -4647,6 +5074,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -4661,6 +5089,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -4676,6 +5105,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) + & *fac_shield(i)*fac_shield(j) enddo return end @@ -7033,6 +7463,252 @@ C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock) return end #endif +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 c------------------------------------------------------------------------------ subroutine eback_sc_corr(esccor) c 7/21/2007 Correlations between the backbone-local and side-chain-local @@ -7174,6 +7850,7 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -7592,6 +8269,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.CONTACTS' include 'COMMON.CHAIN' include 'COMMON.CONTROL' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) integer num_cont_hb_old(maxres) logical lprn,ldone @@ -7894,6 +8572,7 @@ cd write (iout,*) "grij_hb_cont i1",grij_hb_cont(:,jj,i1) call calc_eello(i,jp,i+1,jp1,jj,kk) if (wcorr4.gt.0.0d0) & ecorr=ecorr+eello4(i,jp,i+1,jp1,jj,kk) +CC & *fac_shield(i)**2*fac_shield(j)**2 if (energy_dec.and.wcorr4.gt.0.0d0) 1 write (iout,'(a6,4i5,0pf7.3)') 2 'ecorr4',i,j,i+1,j1,eello4(i,jp,i+1,jp1,jj,kk) @@ -8013,9 +8692,12 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' + include 'COMMON.CONTROL' double precision gx(3),gx1(3) logical lprn lprn=.false. +C print *,"wchodze",fac_shield(i),shield_mode eij=facont_hb(jj,i) ekl=facont_hb(kk,k) ees0pij=ees0p(jj,i) @@ -8024,6 +8706,8 @@ c------------------------------------------------------------------------------ ees0mkl=ees0m(kk,k) ekont=eij*ekl ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl) +C* +C & fac_shield(i)**2*fac_shield(j)**2 cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl) C Following 4 lines for diagnostics. cd ees0pkl=0.0D0 @@ -8036,7 +8720,7 @@ c & ' eij',eij,' eesij',ees0pij,ees0mij,' and ',k,l c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees, c & 'gradcorr_long' C Calculate the multi-body contribution to energy. -c ecorr=ecorr+ekont*ees +C ecorr=ecorr+ekont*ees C Calculate multi-body contributions to the gradient. coeffpees0pij=coeffp*ees0pij coeffmees0mij=coeffm*ees0mij @@ -8087,7 +8771,89 @@ cgrad & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) cgrad enddo cgrad enddo c write (iout,*) "ehbcorr",ekont*ees +C print *,ekont,ees,i,k ehbcorr=ekont*ees +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 return end #ifdef MOMENT @@ -8108,9 +8874,9 @@ C--------------------------------------------------------------------------- & auxmat(2,2) iti1 = itortyp(itype(i+1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1)) + itj1 = itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -8198,16 +8964,16 @@ cd write (iout,*) "a_chujkl",((a_chuj(iii,jjj,kk,k),iii=1,2),jjj=1,2) if (l.eq.j+1) then C parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itype2loc(itype(i)) else - iti=ntortyp + iti=nloctyp endif - itk1=itortyp(itype(k+1)) - itj=itortyp(itype(j)) + itk1=itype2loc(itype(k+1)) + itj=itype2loc(itype(j)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -8351,17 +9117,17 @@ C End vectors else C Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itype2loc(itype(i)) else - iti=ntortyp + iti=nloctyp endif - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + itk1=itype2loc(itype(k+1)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif C A2 kernel(j-1)T A1T call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), @@ -8699,9 +9465,9 @@ cd endif cd write (iout,*) cd & 'EELLO5: Contacts have occurred for peptide groups',i,j, cd & ' and',k,l - itk=itortyp(itype(k)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + itk=itype2loc(itype(k)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) eello5_1=0.0d0 eello5_2=0.0d0 eello5_3=0.0d0 @@ -8770,7 +9536,7 @@ C Cartesian gradient c goto 1112 c1111 continue C Contribution from graph II - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8851,7 +9617,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV cd1110 continue - call transpose2(EE(1,1,itl),auxmat(1,1)) + call transpose2(EE(1,1,l),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8924,7 +9690,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV 1110 continue - call transpose2(EE(1,1,itj),auxmat(1,1)) + call transpose2(EE(1,1,j),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9221,7 +9987,7 @@ C o o o o C C i i C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - itk=itortyp(itype(k)) + itk=itype2loc(itype(k)) s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i)) s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k)) s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k)) @@ -9513,16 +10279,16 @@ C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. iti=itortyp(itype(i)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) + itk=itype2loc(itype(k)) + itk1=itype2loc(itype(k+1)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) @@ -9531,7 +10297,7 @@ C energy moment and not to the cluster cumulant. s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1)) s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9629,24 +10395,24 @@ C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. cd write (2,*) 'eello_graph4: wturn6',wturn6 - iti=itortyp(itype(i)) - itj=itortyp(itype(j)) + iti=itype2loc(itype(i)) + itj=itype2loc(itype(j)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif - itk=itortyp(itype(k)) + itk=itype2loc(itype(k)) if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) + itk1=itype2loc(itype(k+1)) else - itk1=ntortyp + itk1=nloctyp endif - itl=itortyp(itype(l)) + itl=itype2loc(itype(l)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, @@ -9866,11 +10632,11 @@ c j=i+4 k=i+1 l=i+3 - iti=itortyp(itype(i)) - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + iti=itype2loc(itype(i)) + itk=itype2loc(itype(k)) + itk1=itype2loc(itype(k+1)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) cd write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj cd write (2,*) 'i',i,' k',k,' j',j,' l',l cd if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then @@ -10530,11 +11296,12 @@ C first for shielding is setting of function of side-chains 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/ + &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) + &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 @@ -10543,14 +11310,17 @@ 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(k+nres,j)-(c(i,j)+c(i+1,j))/2.0d0 + 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(k+nres,j)-c(k,j) - pept_group(j)=c(i,j)-c(i+1,j) + 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 @@ -10558,8 +11328,14 @@ C lets have their lenght 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 @@ -10567,7 +11343,7 @@ 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)=k + 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 @@ -10577,20 +11353,29 @@ C Lets have the sscale value else scale_fac_dist=-sh_frac_dist*sh_frac_dist & *(2.0*sh_frac_dist-3.0d0) - fac_help_scale=6.0*(scale_fac_dist-scale_fac_dist**2) + 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+short**2/dist_pep_side**2) + costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2) C now costhet_grad - costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side +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 @@ -10602,24 +11387,280 @@ C pep_side0pept_group is vector multiplication do j=1,3 pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j) enddo - fac_alfa_sin=1.0-(pep_side0pept_group/ - & (dist_pep_side*dist_side_calf))**2 + 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 - cosphi=1.0d0/dsqrt(1+rkprim**2/dist_pep_side**2) +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 -C if ((cosphi.le.0.0).or.(costhet.le.0.0)) write(iout,*) "ERROR", -C & cosphi,costhet -C now should be fac_side_grad(k) which will be gradient of factor k which also -C affect the gradient of peptide group i fac_pept_grad(i) and i+1 - write(2,*) "myvolume",VofOverlap,VSolvSphere_div,VolumeTotal - enddo -C write(2,*) "TOTAL VOLUME",i,VolumeTotal -C the scaling factor of the shielding effect + enddo fac_shield(i)=VolumeTotal*div77_81+div4_81 - write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i) +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 return end