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=e86bb6e5de8c25d35efb9c03e5eb07ea42b73499;hb=c1827efa66e69c93c6e2b2e7420b06b430c3550a;hp=6df357a7656a9243d605896014a5a1a41d3c1a70;hpb=7f0180cec09003eb00ab68b95110de7bd15d81a1;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 6df357a..e86bb6e 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -24,6 +24,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.MD' include 'COMMON.CONTROL' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -54,6 +55,8 @@ C FG slaves as WEIGHTS array. weights_(17)=wbond weights_(18)=scal14 weights_(21)=wsccor + weights_(22)=wtube + C FG Master broadcasts the WEIGHTS_ array call MPI_Bcast(weights_(1),n_ene, & MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) @@ -80,6 +83,7 @@ C FG slaves receive the WEIGHTS array wbond=weights(17) scal14=weights(18) wsccor=weights(21) + wtube=weights(22) endif time_Bcast=time_Bcast+MPI_Wtime()-time00 time_Bcastw=time_Bcastw+MPI_Wtime()-time00 @@ -98,6 +102,7 @@ c endif C C Compute the side-chain and electrostatic interaction energy C +C print *,ipot goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj(evdw) @@ -111,6 +116,7 @@ C Berne-Pechukas potential (dilated LJ, angular dependence). goto 107 C Gay-Berne potential (shifted LJ, angular dependence). 104 call egb(evdw) +C print *,"bylem w egb" goto 107 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). 105 call egbv(evdw) @@ -121,6 +127,11 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +cmc +cmc Sep-06: egb takes care of dynamic ss bonds too +cmc +c if (dyn_ss) call dyn_set_nss + c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time01=MPI_Wtime() @@ -129,6 +140,16 @@ c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif +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 c print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then #ifdef SPLITELE @@ -151,9 +172,9 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" eello_turn4=0.0d0 endif else -c write (iout,*) "Soft-spheer ELEC potential" - call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, - & eello_turn4) + write (iout,*) "Soft-spheer ELEC potential" +c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, +c & eello_turn4) endif c print *,"Processor",myrank," computed UELEC" C @@ -185,22 +206,39 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) + 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 endif c print *,"Processor",myrank," computed UB" C C Calculate the SC local energy. C +C print *,"TU DOCHODZE?" call esc(escloc) c print *,"Processor",myrank," computed USC" 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 @@ -209,7 +247,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 @@ -223,6 +261,7 @@ C else esccor=0.0d0 endif +C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C C 12/1/95 Multi-body terms @@ -255,6 +294,26 @@ C after the equilibration time Uconst=0.0d0 Uconst_back=0.0d0 endif +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 + if (TUBElog.gt.0) then +C print *,"just before call" + call calctube(Etube) + else + Etube=0.0d0 + endif + #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -296,8 +355,15 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor + energia(22)=eliptran + energia(23)=Eafmforce + energia(24)=ethetacnstr + energia(25)=Etube +c Here are the energies showed per procesor if the are more processors +c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) + if (dyn_ss) call dyn_set_nss c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -384,20 +450,27 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) + eliptran=energia(22) + Eafmforce=energia(23) + ethetacnstr=energia(24) + Etube=energia(25) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce + & +ethetacnstr+wtube*Etube #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 + & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + & +Eafmforce + & +ethetacnstr+wtube*Etube #endif energia(0)=etot c detecting NaNQ @@ -433,9 +506,10 @@ cMS$ATTRIBUTES C :: proc_proc #endif #ifdef MPI include 'mpif.h' - double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) #endif + double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres), + & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres) + & ,gloc_scbuf(3,-1:maxres) include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.FFIELD' @@ -488,7 +562,7 @@ c enddo call flush(iout) #endif #ifdef SPLITELE - do i=1,nct + do i=0,nct do j=1,3 gradbufc(j,i)=wsc*gvdwc(j,i)+ & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ @@ -499,10 +573,21 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & 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) + & +wtube*gg_tube(j,i) + + + enddo enddo #else - do i=1,nct + do i=0,nct do j=1,3 gradbufc(j,i)=wsc*gvdwc(j,i)+ & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ @@ -514,6 +599,16 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & 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) + & +wtube*gg_tube(j,i) + + + enddo enddo #endif @@ -527,7 +622,7 @@ c enddo enddo call flush(iout) #endif - do i=1,nres + do i=0,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) enddo @@ -570,7 +665,7 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,nnt,-1 + do i=nres-2,-1,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -591,7 +686,7 @@ c enddo enddo call flush(iout) #endif - do i=1,nres + do i=-1,nres do j=1,3 gradbufc_sum(j,i)=gradbufc(j,i) gradbufc(j,i)=0.0d0 @@ -600,7 +695,7 @@ c enddo do j=1,3 gradbufc(j,nres-1)=gradbufc_sum(j,nres) enddo - do i=nres-2,nnt,-1 + do i=nres-2,-1,-1 do j=1,3 gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) enddo @@ -628,9 +723,16 @@ c enddo do k=1,3 gradbufc(k,nres)=0.0d0 enddo - do i=1,nct + 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)+ @@ -649,11 +751,25 @@ c enddo & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & +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) + & +wtube*gg_tube(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)+ @@ -668,12 +784,37 @@ c enddo & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & +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) + & +wtube*gg_tube(j,i) + + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & 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) + & +wtube*gg_tube_sc(j,i) + + + enddo enddo #ifdef DEBUG @@ -708,7 +849,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -717,7 +858,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG do i=1,nres do j=1,3 gloc_scbuf(j,i)=gloc_sc(j,i,icg) @@ -737,7 +878,7 @@ c enddo call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 -#define DEBUG +c#define DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -746,7 +887,7 @@ c enddo enddo enddo #endif -#undef DEBUG +c#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -876,6 +1017,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 @@ -911,6 +1053,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 @@ -962,15 +1109,20 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) + eliptran=energia(22) + Eafmforce=energia(23) + ethetacnstr=energia(24) + etube=energia(25) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, & 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,wsccor, - & edihcnstr,ebr*nss, - & Uconst,etot + & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr, + & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etube,wtube, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -992,9 +1144,14 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ + & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & 'ETOT= ',1pE16.6,' (total)') + #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, & estr,wbond,ebe,wang, @@ -1002,7 +1159,9 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, - & ebr*nss,Uconst,etot + & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etube,wtube, + & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1023,8 +1182,12 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ + & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1079,13 +1242,14 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C have you changed here? + e1=fac*fac*aa + e2=fac*bb evdwij=e1+e2 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') -cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) evdw=evdw+evdwij @@ -1229,8 +1393,9 @@ C rij=1.0D0/r_inv_ij r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C have you changed here? + e1=fac*fac*aa + e2=fac*bb evdwij=e_augm+e1+e2 cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -1356,17 +1521,18 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions C to its derivatives +C have you changed here? fac=(rrij*sigsq)**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa cd write (iout,'(2(a3,i3,2x),15(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & epsi,sigm,chi1,chi2,chip1,chip2, @@ -1410,14 +1576,23 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' logical lprn + integer xshift,yshift,zshift + evdw=0.0D0 ccccc energy_dec=.false. -c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon +C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 +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 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle @@ -1425,6 +1600,70 @@ c if (icall.eq.0) lprn=.false. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+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 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 + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1437,6 +1676,38 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + 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 ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1465,17 +1736,119 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi + 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 write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) +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 dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) +C xj=xj-xi +C yj=yj-yi +C zj=zj-zi c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi c write (iout,*) "j",j," dc_norm", c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + 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 C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -1496,18 +1869,24 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C here to start with +C if (c(i,3).gt. + faclip=fac + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt +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) c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij + evdw=evdw+evdwij*sss if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,chi1,chi2,chip1,chip2, @@ -1524,16 +1903,32 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac +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 c fac=0.0d0 C Calculate the radial part of the gradient + 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 gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad + endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -1570,6 +1965,41 @@ c if (icall.eq.0) lprn=.true. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) + 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 + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1606,9 +2036,75 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi +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)) +C write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 + xj=xj_safe+xshift*boxxsize + yj=yj_safe+yshift*boxysize + zj=zj_safe+zshift*boxzsize + dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -1629,8 +2125,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa + e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt @@ -1639,8 +2135,8 @@ c--------------------------------------------------------------- evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij+e_augm if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi),i,restyp(itypj),j, & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0), @@ -1654,6 +2150,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac-2*expon*rrij*e_augm + fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -1741,6 +2238,7 @@ C---------------------------------------------------------------------------- include 'COMMON.CALC' include 'COMMON.IOUNITS' double precision dcosom1(3),dcosom2(3) +cc print *,'sss=',sss eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 @@ -1759,16 +2257,16 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k)) enddo do k=1,3 - gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k) + gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss enddo c write (iout,*) "gg",(gg(k),k=1,3) do k=1,3 - gvdwx(k,i)=gvdwx(k,i)-gg(k) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k) & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) - & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv - gvdwx(k,j)=gvdwx(k,j)+gg(k) + & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k) & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) - & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) @@ -1783,8 +2281,8 @@ cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad enddo do l=1,3 - gvdwc(l,i)=gvdwc(l,i)-gg(l) - gvdwc(l,j)=gvdwc(l,j)+gg(l) + gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l) enddo return end @@ -1886,7 +2384,7 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' dimension ggg(3) -cd write(iout,*) 'In EELEC_soft_sphere' +C write(iout,*) 'In EELEC_soft_sphere' ees=0.0D0 evdw1=0.0D0 eel_loc=0.0d0 @@ -1901,6 +2399,12 @@ cd write(iout,*) 'In EELEC_soft_sphere' xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -1914,10 +2418,49 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) dxj=dc(1,j) dyj=dc(2,j) dzj=dc(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - 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 + 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 rij=xj*xj+yj*yj+zj*zj + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) if (rij.lt.r0ijsq) then evdw1ij=0.25d0*(rij-r0ijsq)**2 fac=rij-r0ijsq @@ -1925,13 +2468,13 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) evdw1ij=0.0d0 fac=0.0d0 endif - evdw1=evdw1+evdw1ij + evdw1=evdw1+evdw1ij*sss C C Calculate contributions to the Cartesian gradient. C - ggg(1)=fac*xj - ggg(2)=fac*yj - ggg(3)=fac*zj + ggg(1)=fac*xj*sssgrad + ggg(2)=fac*yj*sssgrad + ggg(3)=fac*zj*sssgrad do k=1,3 gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) @@ -2250,6 +2793,98 @@ C C Compute the virtual-bond-torsional-angle dependent quantities needed C to calculate the el-loc multibody terms of various order. C +c write(iout,*) 'nphi=',nphi,nres +#ifdef PARMAT + do i=ivec_start+2,ivec_end+2 +#else + do i=3,nres+1 +#endif +#ifdef NEWCORR + 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 +c write(iout,*),i + 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.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.0d0) + & +bnew2(2,1,iti)*dsin(theta(i-1)) + & +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 + & +bnew2(2,1,iti)*dcos(theta(i-1)) + & -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 +c if (ggb1(1,i).eq.0.0d0) then +c write(iout,*) 'i=',i,ggb1(1,i), +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 endif + b1(2,i-2)=bnew1(1,2,iti) + gtb1(2,i-2)=0.0 + b2(2,i-2)=bnew2(1,2,iti) + gtb2(2,i-2)=0.0 + EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1)) + 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) + gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1)) + gtEE(1,2,i-2)=0.0d0 + gtEE(2,2,i-2)=0.0d0 + gtEE(2,1,i-2)=0.0d0 +c EE(2,2,iti)=0.0d0 +c EE(1,2,iti)=0.5d0*eenew(1,iti) +c EE(2,1,iti)=0.5d0*eenew(1,iti) +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)) + 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) +c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) +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) + 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 @@ -2323,15 +2958,15 @@ C 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+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 cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) @@ -2339,8 +2974,18 @@ cd write (iout,*) 'b2',b2(:,iti) cd write (iout,*) 'Ug',Ug(:,:,i-2) c if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then - call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2)) - call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2)) + call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) +#ifdef NEWCORR + 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,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", +c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2), +c & eug(2,2,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2)) @@ -2362,21 +3007,32 @@ c if (i .gt. iatel_s+2) then enddo enddo endif - call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2)) - call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2)) + call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2)) + call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2)) do k=1,2 muder(k,i-2)=Ub2der(k,i-2) enddo 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)) + if (itype(i-1).le.ntyp) then + iti1 = itype2loc(itype(i-1)) + else + iti1=nloctyp + endif else - iti1=ntortyp+1 + iti1=nloctyp endif do k=1,2 - mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) + mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -cd write (iout,*) 'mu ',mu(:,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) @@ -2387,7 +3043,7 @@ cd write (iout,*) 'mu2',mu2(:,i-2) call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2)) call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2)) C Vectors and matrices dependent on a single virtual-bond dihedral. - call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) + call matvec2(DD(1,1,iti),b1tilde(1,i-1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2)) @@ -2663,11 +3319,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 @@ -2700,10 +3356,11 @@ C include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -2782,9 +3439,25 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms C 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 +c if (i.le.1) cycle +C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle +C changes suggested by Ana to avoid out of bounds +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 +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) @@ -2794,15 +3467,30 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + xmedi=mod(xmedi,boxxsize) + if (xmedi.lt.0) xmedi=xmedi+boxxsize + ymedi=mod(ymedi,boxysize) + if (ymedi.lt.0) ymedi=ymedi+boxysize + zmedi=mod(zmedi,boxzsize) + if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end + if (i.lt.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 +C changes suggested by Ana to avoid out of bounds +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) cycle + & .or. itype(i+4).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) dzi=dc(3,i) @@ -2812,17 +3500,64 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi +C Return atom into box, boxxsize is size of box in x dimension +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 + num_conti=num_cont_hb(i) +c write(iout,*) "JESTEM W PETLI" call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i +C Loop over all neighbouring boxes +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c +CTU KURWA do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) 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 +c & .or.((i+2).gt.nres) +c & .or.((i-1).le.0) +C end of changes by Ana +c & .or. itype(i+2).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2832,15 +3567,65 @@ c xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi + 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 + 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 write (iout,*) i,j,itype(i),itype(j) - if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle +C do j=16,17 +C write (iout,*) i,j +C 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 +c & .or.((j+2).gt.nres) +c & .or.((j-1).le.0) +C end of changes by Ana +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 num_cont_hb(i)=num_conti enddo ! i +C enddo ! zshift +C enddo ! yshift +C enddo ! xshift + c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres cd write (iout,'(i3,3f10.5,5x,3f10.5)') @@ -2871,10 +3656,13 @@ C------------------------------------------------------------------------------- include 'COMMON.VECTORS' include 'COMMON.FFIELD' include 'COMMON.TIME1' + include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), - & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) + & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), + & gmuij2(4),gmuji2(4) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 @@ -2889,6 +3677,7 @@ C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ + integer xshift,yshift,zshift c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j c ind=ind+1 @@ -2905,10 +3694,84 @@ c ind=ind+1 dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi +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 rij=xj*xj+yj*yj+zj*zj + + sss=sscale(sqrt(rij)) + sssgrad=sscagrad(sqrt(rij)) +c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij @@ -2924,50 +3787,146 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions ev2=bbb*r6ij fac3=ael6i*r6ij fac4=ael3i*r3ij - evdwij=ev1+ev2 + evdwij=(ev1+ev2) el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)) el2=fac4*fac - eesij=el1+el2 +C MARYSIA +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 +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 - evdw1=evdw1+evdwij + else + fac_shield(i)=1.0 + fac_shield(j)=1.0 + eesij=(el1+el2) + ees=ees+eesij + endif + evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij - write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij + write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') + &'evdw1',i,j,evdwij + &,iteli,itelj,aaa,evdw1 + write (iout,*) sss + write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, + &fac_shield(i),fac_shield(j) endif C C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij) + facvdw=-6*rrmij*(ev1+evdwij)*sss facel=-3*rrmij*(el1+eesij) fac1=fac 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. * @@ -2976,9 +3935,15 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + if (sss.gt.0.0) then + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + else + ggg(1)=0.0 + ggg(2)=0.0 + ggg(3)=0.0 + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -2998,8 +3963,9 @@ cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) cgrad enddo cgrad enddo #else - facvdw=ev1+evdwij - facel=el1+eesij +C MARYSIA + facvdw=(ev1+evdwij)*sss + facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij @@ -3009,8 +3975,11 @@ cgrad enddo * 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 @@ -3030,9 +3999,9 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -3053,7 +4022,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) @@ -3069,16 +4039,23 @@ 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 & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN @@ -3089,6 +4066,7 @@ C Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al. C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms C are computed for EVERY pair of non-contiguous peptide groups. C + if (j.lt.nres-1) then j1=j+1 j2=j-1 @@ -3097,10 +4075,20 @@ C j2=j-2 endif kkk=0 + lll=0 do k=1,2 do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) +c write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l +#ifdef NEWCORR + gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) +c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) + gmuij2(kkk)=gUb2(k,i)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) +c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) + gmuji2(kkk)=mu(k,i)*gUb2(l,j) +#endif enddo enddo cd write (iout,*) 'EELEC: i',i,' j',j @@ -3266,24 +4254,138 @@ 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=',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) + & +a23*gmuij1(2) + & +a32*gmuij1(3) + & +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) + gloc(nphi+i,icg)=gloc(nphi+i,icg)+ + & geel_loc_ij*wel_loc +c write(iout,*) "derivative over thatai-1" +c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), +c & a33*gmuij2(4) + geel_loc_ij= + & a22*gmuij2(1) + & +a23*gmuij2(2) + & +a32*gmuij2(3) + & +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) + & +a32*gmuji1(3) + & +a33*gmuji1(4) +c write(iout,*) "derivative over thataj" +c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), +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) + & +a32*gmuji2(3) + & +a33*gmuji2(4) +c write(iout,*) "derivative over thataj-1" +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 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij +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) eel_loc=eel_loc+eel_loc_ij 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) @@ -3297,14 +4399,22 @@ cgrad enddo cgrad enddo 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) - 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) + 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 @@ -3383,8 +4493,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 @@ -3452,17 +4572,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 @@ -3514,10 +4646,13 @@ 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), - & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2) + & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2), + & gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2), + & auxgmat2(2,2),auxgmatt2(2,2) double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, @@ -3541,11 +4676,83 @@ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn3(i,a_temp,eello_turn3_num) call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) +c auxalary matices for theta gradient +c auxalary matrix for i+1 and constant i+2 + call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1)) +c auxalary matrix for i+2 and constant i+1 + call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) + call transpose2(auxgmat1(1,1),auxgmatt1(1,1)) + call transpose2(auxgmat2(1,1),auxgmatt2(1,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)) - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'eturn3',i,j,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#ifdef NEWCORR +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) +C#endif + +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 @@ -3554,12 +4761,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) @@ -3573,6 +4782,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) @@ -3580,6 +4791,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 @@ -3587,6 +4799,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) @@ -3594,6 +4807,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 @@ -3614,10 +4828,15 @@ 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), - & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2) + & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2), + & auxgEvec1(2),auxgEvec2(2),auxgEvec3(2), + & gte1t(2,2),gte2t(2,2),gte3t(2,2), + & gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2), + & gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2) double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, @@ -3637,58 +4856,192 @@ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn4(i,a_temp,eello_turn4_num) c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 +c write(iout,*)"WCHODZE W PROGRAM" a_temp(1,1)=a22 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)) call transpose2(Eug(1,1,i+3),e3t(1,1)) +C Ematrix derivative in theta + call transpose2(gtEUg(1,1,i+1),gte1t(1,1)) + call transpose2(gtEug(1,1,i+2),gte2t(1,1)) + call transpose2(gtEug(1,1,i+3),gte3t(1,1)) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) +c eta1 in derivative theta + call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) +c auxgvec is derivative of Ub2 so i+3 theta + call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) +c auxalary matrix of E i+1 + call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1)) +c s1=0.0 +c gs1=0.0 + s1=scalar2(b1(1,i+2),auxvec(1)) +c derivative of theta i+2 with constant i+3 + gs23=scalar2(gtb1(1,i+2),auxvec(1)) +c derivative of theta i+2 with constant i+2 + gs32=scalar2(b1(1,i+2),auxgvec(1)) +c derivative of E matix in theta of i+1 + gsE13=scalar2(b1(1,i+2),auxgEvec1(1)) + call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) +c ea31 in derivative theta + call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) +c auxilary matrix auxgvec of Ub2 with constant E matirx + call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1)) +c auxilary matrix auxgEvec1 of E matix with Ub2 constant + call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1)) + +c s2=0.0 +c gs2=0.0 + s2=scalar2(b1(1,i+1),auxvec(1)) +c derivative of theta i+1 with constant i+3 + gs13=scalar2(gtb1(1,i+1),auxvec(1)) +c derivative of theta i+2 with constant i+1 + gs21=scalar2(b1(1,i+1),auxgvec(1)) +c derivative of theta i+3 with constant i+1 + gsE31=scalar2(b1(1,i+1),auxgEvec3(1)) +c write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2), +c & gtb1(1,i+1) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) +c two derivatives over diffetent matrices +c gtae3e2 is derivative over i+3 + call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1)) +c ae3gte2 is derivative over i+2 + call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) +c three possible derivative over theta E matices +c i+1 + call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1)) +c i+2 + call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1)) +c i+3 + call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) + + 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) - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'eturn4',i,j,-(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 + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') + & 'eturn4',i,j,-(s1+s2+s3) +c write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), +c & ' eello_turn4_num',8*eello_turn4_num C Derivatives in gamma(i) call transpose2(EUgder(1,1,i+1),e1tder(1,1)) call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1)) call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) 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)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,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)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1)) call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1)) 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 @@ -3699,15 +5052,16 @@ C Derivatives of this turn contributions in DC(i+2) a_temp(2,2)=agg(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) 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 @@ -3718,57 +5072,61 @@ C Remaining derivatives of this turn contribution a_temp(2,2)=aggi(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) 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) a_temp(2,2)=aggi1(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) 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) a_temp(2,2)=aggj(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) 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) a_temp(2,2)=aggj1(l,4) call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) - s1=scalar2(b1(1,iti2),auxvec(1)) + s1=scalar2(b1(1,i+2),auxvec(1)) call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) - s2=scalar2(b1(1,iti1),auxvec(1)) + s2=scalar2(b1(1,i+1),auxvec(1)) call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) 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 @@ -3829,13 +5187,49 @@ C r0_scp=4.5d0 cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - +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 do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -3846,10 +5240,75 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + 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 rij=xj*xj+yj*yj+zj*zj + r0ij=r0_scp r0ijsq=r0ij*r0ij if (rij.lt.r0ijsq) then @@ -3900,6 +5359,9 @@ cgrad enddo enddo ! iint enddo ! i +C enddo !zshift +C enddo !yshift +C enddo !xshift return end C----------------------------------------------------------------------------- @@ -3920,18 +5382,59 @@ C include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' dimension ggg(3) evdw2=0.0D0 evdw2_14=0.0d0 +c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e +C do xshift=-1,1 +C do yshift=-1,1 +C do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize +c xi=xi+xshift*boxxsize +c yi=yi+yshift*boxysize +c zi=zi+zshift*boxzsize +c print *,xi,yi,zi,'polozenie i' +C Return atom into box, boxxsize is size of box in x dimension +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" +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 do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -3942,26 +5445,97 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + 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' rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c print *,rrij + sss=sscale(1.0d0/(dsqrt(rrij))) +c print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz' +c if (sss.eq.0) print *,'czasem jest OK' + if (sss.le.0.0d0) cycle + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) e2=fac*bad(itypj,iteli) if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+e1+e2 + evdw2_14=evdw2_14+(e1+e2)*sss endif evdwij=e1+e2 - evdw2=evdw2+evdwij - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw2',i,j,evdwij + evdw2=evdw2+evdwij*sss + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') + & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), + & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C - fac=-(evdwij+e1)*rrij + fac=-(evdwij+e1)*rrij*sss + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -3996,10 +5570,14 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - enddo +c endif !endif for sscale cutoff + enddo ! j enddo ! iint enddo ! i +c enddo !zshift +c enddo !yshift +c enddo !xshift do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) @@ -4031,8 +5609,13 @@ C include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 + do i=1,3 + ggg(i)=0.0d0 + enddo +C write (iout,*) ,"link_end",link_end,constr_dist cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr cd write(iout,*)'link_start=',link_start,' link_end=',link_end if (link_end.eq.0) return @@ -4049,53 +5632,118 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +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. & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij +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 else C Calculate the distance between the two points and its difference from the C target distance. - dd=dist(ii,jj) - rdis=dd-dhpb(i) + 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) C Get the force constant corresponding to this distance. - waga=forcon(i) + waga=forcon(i) C Calculate the contribution to energy. - ehpb=ehpb+waga*rdis*rdis + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C - 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 + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) C If this is a SC-SC distance, we need to calculate the contributions to the C Cartesian gradient in the SC vectors (ghpbx). - if (iii.lt.ii) then + if (iii.lt.ii) then do j=1,3 ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo - endif + endif cgrad do j=iii,jjj-1 cgrad do k=1,3 cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) cgrad enddo cgrad enddo - 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 endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -4155,7 +5803,7 @@ c dscj_inv=dsc_inv(itypj) cosphi=om12-om1*om2 eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) & +akct*deltad*deltat12 - & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi + & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, c & " deltat12",deltat12," eij",eij @@ -4208,25 +5856,34 @@ c estr=0.0d0 estr1=0.0d0 do i=ibondp_start,ibondp_end - 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 + 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 + if (energy_dec) write(iout,*) "dum_bond",i,diff + else +C NO vbldp0 is the equlibrium lenght of spring for peptide group diff = vbld(i)-vbldp0 - if (energy_dec) write (iout,*) + endif + if (energy_dec) write (iout,'(a7,i5,4f7.3)') & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) enddo c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) - endif +c endif enddo + estr=0.5d0*AKP*estr+estr1 c c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included @@ -4237,7 +5894,7 @@ c nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) - if (energy_dec) write (iout,*) + if (energy_dec) write (iout,*) & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, & AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff @@ -4279,7 +5936,7 @@ c end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -4296,6 +5953,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@ -4306,7 +5964,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) @@ -4323,7 +5982,7 @@ C Zero the energy function and its derivative at 0 or pi. ichir22=isign(1,itype(i)) endif - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4336,7 +5995,7 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4344,8 +6003,8 @@ C Zero the energy function and its derivative at 0 or pi. z(1)=cos(phii1) #else phii1=phi(i+1) - z(1)=dcos(phii1) #endif + z(1)=dcos(phii1) z(2)=dsin(phii1) else z(1)=0.0D0 @@ -4363,6 +6022,7 @@ C In following comments this theta will be referred to as t_c. bthetk=bthet(k,itype2,ichir21,ichir22) endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) +c write(iout,*) 'chuj tu', y(k),z(k) enddo dthett=thet_pred_mean*ssd thet_pred_mean=thet_pred_mean*ss+a0thet(it) @@ -4399,12 +6059,40 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai - if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'ebend',i,ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)') + & 'ebend',i,ethetai,theta(i),itype(i) if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 - gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg) enddo + 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 + C Ufff.... We've done all this!!! return end @@ -4421,7 +6109,8 @@ C--------------------------------------------------------------------------- C Calculate the contributions to both Gaussian lobes. C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time) C The "polynomial part" of the "standard deviation" of this part of -C the distribution. +C the distributioni. +ccc write (iout,*) thetai,thet_pred_mean sig=polthet(3,it) do j=2,0,-1 sig=sig*thet_pred_mean+polthet(j,it) @@ -4451,6 +6140,7 @@ C Following variable is sigma(t_c)**(-2) delthe0=thetai-theta0i term1=-0.5D0*sigcsq*delthec*delthec term2=-0.5D0*sig0inv*delthe0*delthe0 +C write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0 C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and C NaNs in taking the logarithm. We extract the largest exponent which is added C to the energy (this being the log of the distribution) at the end of energy @@ -4478,6 +6168,7 @@ C Contribution of the bending energy from this theta is just the -log of C the sum of the contributions from the two lobes and the pre-exponential C factor. Simple enough, isn't it? ethetai=(-dlog(termexp)-termm+dlog(termpre)) +C write (iout,*) 'termexp',termexp,termm,termpre,i C NOW the derivatives!!! C 6/6/97 Take into account the deformation. E_theta=(delthec*sigcsq*term1 @@ -4518,7 +6209,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -4537,6 +6228,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -4544,7 +6236,10 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle +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 (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -4556,7 +6251,8 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then +C print *,ethetai + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4571,13 +6267,13 @@ C propagation of chirality for glycine type enddo else phii=0.0d0 - ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4592,7 +6288,7 @@ C propagation of chirality for glycine type enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -4611,14 +6307,14 @@ C propagation of chirality for glycine type sinph1ph2(k,l)=scl-csl enddo enddo -c if (lprn) then + if (lprn) then write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2, & " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1 -c write (iout,*) "coskt and sinkt" -c do k=1,nntheterm -c write (iout,*) k,coskt(k),sinkt(k) -c enddo -c endif + write (iout,*) "coskt and sinkt" + do k=1,nntheterm + write (iout,*) k,coskt(k),sinkt(k) + enddo + endif do k=1,ntheterm ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k) dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock) @@ -4642,6 +6338,7 @@ c endif enddo write(iout,*) "ethetai",ethetai endif +C print *,ethetai do m=1,ntheterm2 do k=1,nsingle aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k) @@ -4662,10 +6359,16 @@ c endif & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet", & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet", & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) enddo enddo +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) if (lprn) & write(iout,*) "ethetai",ethetai +C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) do m=1,ntheterm3 do k=2,ndouble do l=1,k-1 @@ -4701,16 +6404,46 @@ c endif enddo 10 continue c lprn1=.true. -c if (lprn1) - write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') +C print *,ethetai + if (lprn1) + & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai c lprn1=.false. etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo +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 + return end #endif @@ -5048,7 +6781,7 @@ C cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) + it=iabs(itype(i)) if (it.eq.10) goto 1 c C Compute the axes of tghe local cartesian coordinates system; store in @@ -5066,7 +6799,7 @@ C & dc_norm(3,i+nres) y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac enddo do j = 1,3 - z_prime(j) = -uz(j,i-1)*dsign(1.0,dfloat(itype(i))) + z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i))) enddo c write (2,*) "i",i c write (2,*) "x_prime",(x_prime(j),j=1,3) @@ -5106,7 +6839,7 @@ c write (2,*) "xx",xx," yy",yy," zz",zz Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@ -5148,7 +6881,8 @@ c & sumene4, c & dscp1,dscp2,sumene c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene -c write (2,*) "i",i," escloc",sumene,escloc +c write (2,*) "i",i," escloc",sumene,escloc,it,itype(i) +c & ,zz,xx,yy c#define DEBUG #ifdef DEBUG C @@ -5280,9 +7014,9 @@ c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i) dZZ_Ci(k)=0.0d0 do j=1,3 dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) - & *dsign(1.0,dfloat(itype(i)))*dC_norm(j,i+nres) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) - & *dsign(1.0,dfloat(itype(i)))*dC_norm(j,i+nres) + & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) enddo dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres)) @@ -5474,9 +7208,9 @@ c lprn=.true. do i=iphi_start,iphi_end etors_ii=0.0D0 if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Proline-Proline pair is a special case... @@ -5529,12 +7263,12 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + 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*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@ -5570,8 +7304,15 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle +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 etors_ii=0.0D0 if (iabs(itype(i)).eq.20) then iblock=2 @@ -5633,18 +7374,21 @@ c do i=1,ndih_constr difi=pinorm(phii-phi0(i)) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + 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*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + 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 -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) + 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 enddo cd write (iout,*) 'edihcnstr',edihcnstr return @@ -5672,8 +7416,15 @@ c lprn=.true. etors_d=0.0D0 c write(iout,*) "a tu??" do i=iphid_start,iphid_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle +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 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5683,9 +7434,21 @@ c write(iout,*) "a tu??" gloci2=0.0D0 iblock=1 if (iabs(itype(i+1)).eq.20) iblock=2 +C Iblock=2 Proline type +C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT +C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO- +C if (itype(i+1).eq.ntyp1) iblock=3 +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 C Regular cosine and sine terms do j=1,ntermd_1(itori,itori1,itori2,iblock) +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) v1cij=v1c(1,j,itori,itori1,itori2,iblock) v1sij=v1s(1,j,itori,itori1,itori2,iblock) v2cij=v1c(2,j,itori,itori1,itori2,iblock) @@ -5723,6 +7486,252 @@ C Regular cosine and sine terms 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 @@ -5864,6 +7873,7 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -6282,6 +8292,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 @@ -6584,6 +8595,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) @@ -6703,9 +8715,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) @@ -6714,6 +8729,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 @@ -6726,7 +8743,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 @@ -6777,7 +8794,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 @@ -6798,17 +8897,17 @@ 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+1 + itj1=nloctyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) dipderi(iii)=Ub2der(iii,i) - dipi(iii,2)=b1(iii,iti1) + dipi(iii,2)=b1(iii,i+1) dipj(iii,1)=Ub2(iii,j) dipderj(iii)=Ub2der(iii,j) - dipj(iii,2)=b1(iii,itj1) + dipj(iii,2)=b1(iii,j+1) enddo kkk=0 do iii=1,2 @@ -6888,16 +8987,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+1 + 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+1 + itl1=nloctyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -6988,26 +9087,26 @@ C They are needed only when the fifth- or the sixth-order cumulants are C indluded. IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN call transpose2(AEA(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1)) call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1)) call transpose2(AEAderg(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1)) - call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1)) - call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1)) + call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1)) + call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1)) call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1)) call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1)) call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1)) call transpose2(AEA(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),b1(1,j),AEAb1(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2)) call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2)) call transpose2(AEAderg(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2)) + call matvec2(auxmat(1,1),b1(1,j),AEAb1derg(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2)) - call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2)) - call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2)) + call matvec2(AEA(1,1,2),b1(1,l+1),AEAb1(1,2,2)) + call matvec2(AEAderg(1,1,2),b1(1,l+1),AEAb1derg(1,2,2)) call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2)) call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2)) call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2)) @@ -7016,20 +9115,20 @@ C Calculate the Cartesian derivatives of the vectors. do kkk=1,5 do lll=1,3 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti), + call matvec2(auxmat(1,1),b1(1,i), & AEAb1derx(1,lll,kkk,iii,1,1)) call matvec2(auxmat(1,1),Ub2(1,i), & AEAb2derx(1,lll,kkk,iii,1,1)) - call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1), & AEAb1derx(1,lll,kkk,iii,2,1)) call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1), & AEAb2derx(1,lll,kkk,iii,2,1)) call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj), + call matvec2(auxmat(1,1),b1(1,j), & AEAb1derx(1,lll,kkk,iii,1,2)) call matvec2(auxmat(1,1),Ub2(1,j), & AEAb2derx(1,lll,kkk,iii,1,2)) - call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1), + call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,l+1), & AEAb1derx(1,lll,kkk,iii,2,2)) call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1), & AEAb2derx(1,lll,kkk,iii,2,2)) @@ -7041,17 +9140,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+1 + 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+1 + 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), @@ -7126,26 +9225,26 @@ C indluded. IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or. & (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN call transpose2(AEA(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1)) call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1)) call transpose2(AEAderg(1,1,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1)) + call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1)) call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1)) - call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1)) - call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1)) + call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1)) + call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1)) call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1)) call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1)) call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1)) call transpose2(AEA(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),b1(1,j+1),AEAb1(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2)) call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2)) call transpose2(AEAderg(1,1,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2)) + call matvec2(auxmat(1,1),b1(1,l),AEAb1(1,1,2)) call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2)) - call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2)) - call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2)) + call matvec2(AEA(1,1,2),b1(1,j+1),AEAb1(1,2,2)) + call matvec2(AEAderg(1,1,2),b1(1,j+1),AEAb1derg(1,2,2)) call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2)) call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2)) call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2)) @@ -7154,20 +9253,20 @@ C Calculate the Cartesian derivatives of the vectors. do kkk=1,5 do lll=1,3 call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,iti), + call matvec2(auxmat(1,1),b1(1,i), & AEAb1derx(1,lll,kkk,iii,1,1)) call matvec2(auxmat(1,1),Ub2(1,i), & AEAb2derx(1,lll,kkk,iii,1,1)) - call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1), & AEAb1derx(1,lll,kkk,iii,2,1)) call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1), & AEAb2derx(1,lll,kkk,iii,2,1)) call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1)) - call matvec2(auxmat(1,1),b1(1,itl), + call matvec2(auxmat(1,1),b1(1,l), & AEAb1derx(1,lll,kkk,iii,1,2)) call matvec2(auxmat(1,1),Ub2(1,l), & AEAb2derx(1,lll,kkk,iii,1,2)) - call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1), + call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,j+1), & AEAb1derx(1,lll,kkk,iii,2,2)) call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j), & AEAb2derx(1,lll,kkk,iii,2,2)) @@ -7389,9 +9488,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 @@ -7460,11 +9559,11 @@ 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) - eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk)) + eello5_2=scalar2(AEAb1(1,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k)) C Explicit gradient in virtual-dihedral angles. g_corr5_loc(k-1)=g_corr5_loc(k-1) @@ -7474,11 +9573,11 @@ C Explicit gradient in virtual-dihedral angles. vv(2)=pizda(2,1)-pizda(1,2) if (l.eq.j+1) then g_corr5_loc(l-1)=g_corr5_loc(l-1) - & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk)) + & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k))) else g_corr5_loc(j-1)=g_corr5_loc(j-1) - & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk)) + & +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k))) endif C Cartesian gradient @@ -7490,7 +9589,7 @@ C Cartesian gradient vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) derx(lll,kkk,iii)=derx(lll,kkk,iii) - & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk)) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,k)) & -0.5d0*scalar2(vv(1),Ctobr(1,k)) enddo enddo @@ -7541,11 +9640,11 @@ 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) - eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl)) + eello5_4=scalar2(AEAb1(1,2,2),b1(1,l)) & -0.5d0*scalar2(vv(1),Ctobr(1,l)) C Explicit gradient in virtual-dihedral angles. g_corr5_loc(l-1)=g_corr5_loc(l-1) @@ -7554,7 +9653,7 @@ C Explicit gradient in virtual-dihedral angles. vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) g_corr5_loc(k-1)=g_corr5_loc(k-1) - & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl)) + & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,l)) & -0.5d0*scalar2(vv(1),Ctobr(1,l))) C Cartesian gradient do iii=1,2 @@ -7565,7 +9664,7 @@ C Cartesian gradient vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) derx(lll,kkk,iii)=derx(lll,kkk,iii) - & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl)) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,l)) & -0.5d0*scalar2(vv(1),Ctobr(1,l)) enddo enddo @@ -7614,11 +9713,11 @@ 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) - eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj)) + eello5_4=scalar2(AEAb1(1,2,2),b1(1,j)) & -0.5d0*scalar2(vv(1),Ctobr(1,j)) C Explicit gradient in virtual-dihedral angles. g_corr5_loc(j-1)=g_corr5_loc(j-1) @@ -7627,7 +9726,7 @@ C Explicit gradient in virtual-dihedral angles. vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) g_corr5_loc(k-1)=g_corr5_loc(k-1) - & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj)) + & +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,j)) & -0.5d0*scalar2(vv(1),Ctobr(1,j))) C Cartesian gradient do iii=1,2 @@ -7638,7 +9737,7 @@ C Cartesian gradient vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii) - & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj)) + & +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,j)) & -0.5d0*scalar2(vv(1),Ctobr(1,j)) enddo enddo @@ -7702,9 +9801,9 @@ cd ghalf=0.0d0 cold ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k) cgrad ghalf=0.5d0*ggg2(ll) cd ghalf=0.0d0 - gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2) + gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2) gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2) - gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2) + gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2) gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2) gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl @@ -7911,7 +10010,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)) @@ -7920,8 +10019,8 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC vv1(1)=pizda1(1,1)-pizda1(2,2) vv1(2)=pizda1(1,2)+pizda1(2,1) s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i)) - vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk) - vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk) + vv(1)=AEAb1(1,2,imat)*b1(1,k)-AEAb1(2,2,imat)*b1(2,k) + vv(2)=AEAb1(1,2,imat)*b1(2,k)+AEAb1(2,2,imat)*b1(1,k) s5=scalar2(vv(1),Dtobr2(1,i)) cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5) @@ -7934,8 +10033,8 @@ cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1)) vv1(1)=pizda1(1,1)-pizda1(2,2) vv1(2)=pizda1(1,2)+pizda1(2,1) - vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk) - vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk) + vv(1)=AEAb1derg(1,2,imat)*b1(1,k)-AEAb1derg(2,2,imat)*b1(2,k) + vv(2)=AEAb1derg(1,2,imat)*b1(2,k)+AEAb1derg(2,2,imat)*b1(1,k) if (l.eq.j+1) then g_corr6_loc(l-1)=g_corr6_loc(l-1) & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i)) @@ -7974,10 +10073,10 @@ cd write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5 vv1(1)=pizda1(1,1)-pizda1(2,2) vv1(2)=pizda1(1,2)+pizda1(2,1) s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i)) - vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk) - & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk) - vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk) - & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk) + vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,k) + & -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,k) + vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,k) + & +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,k) s5=scalar2(vv(1),Dtobr2(1,i)) derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5) enddo @@ -7999,7 +10098,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -8009,12 +10108,12 @@ C C C o o C C \ /l\ /j\ / C C \ / \ / \ / C -C o| o | | o |o C +C o| o | | o |o C C \ j|/k\| \ |/k\|l C C \ / \ \ / \ C C o o C -C i i C -C C +C i i C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l C AL 7/4/01 s1 would occur in the sixth-order moment, @@ -8185,10 +10284,10 @@ c---------------------------------------------------------------------------- double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C -C o o C +C o o C C /l\ / \ /j\ C C / \ / \ / \ C C /| o |o o| o |\ C @@ -8203,25 +10302,25 @@ 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+1 + 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+1 + itl1=nloctyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) #endif - call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1)) - s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) - call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1)) - s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) - call transpose2(EE(1,1,itk),auxmat(1,1)) + call matvec2(AECA(1,1,1),b1(1,k+1),auxvec(1)) + 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,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) @@ -8235,13 +10334,13 @@ cd & "sum",-(s2+s3+s4) #endif c eello6_graph3=-s4 C Derivatives in gamma(k-1) - call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1)) - s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + call matvec2(AECAderg(1,1,2),b1(1,l+1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k)) g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4) C Derivatives in gamma(l-1) - call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1)) - s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) + call matvec2(AECAderg(1,1,1),b1(1,k+1),auxvec(1)) + s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -8258,12 +10357,12 @@ C Cartesian derivatives. s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k) endif #endif - call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1), + call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,k+1), & auxvec(1)) - s2=0.5d0*scalar2(b1(1,itk),auxvec(1)) - call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1), + s2=0.5d0*scalar2(b1(1,k),auxvec(1)) + call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,l+1), & auxvec(1)) - s3=0.5d0*scalar2(b1(1,itj1),auxvec(1)) + s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1), & pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) @@ -8302,7 +10401,7 @@ c---------------------------------------------------------------------------- & auxvec1(2),auxmat1(2,2) logical swap CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC -C C +C C C Parallel Antiparallel C C C C o o C @@ -8310,33 +10409,33 @@ C /l\ / \ /j\ C C / \ / \ / \ C C /| o |o o| o |\ C C \ j|/k\| \ |/k\|l C -C \ / \ \ / \ C +C \ / \ \ / \ C C o \ o \ C C i i C -C C +C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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+1 + 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+1 + 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+1 + 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, @@ -8351,11 +10450,11 @@ cd & ' itl',itl,' itl1',itl1 call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1)) s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) if (j.eq.l+1) then - call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + call matvec2(ADtEA1(1,1,3-imat),b1(1,j+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec1(1)) else - call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + call matvec2(ADtEA1(1,1,3-imat),b1(1,l+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec1(1)) endif call transpose2(EUg(1,1,k),auxmat(1,1)) call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1)) @@ -8379,11 +10478,11 @@ C Derivatives in gamma(i-1) #endif s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1)) if (j.eq.l+1) then - call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,j+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec1(1)) else - call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,l+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec1(1)) endif s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i)) if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then @@ -8412,11 +10511,11 @@ C Derivatives in gamma(k-1) call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1)) s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1)) if (j.eq.l+1) then - call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,j+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec1(1)) else - call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1)) + call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,l+1),auxvec1(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec1(1)) endif call transpose2(EUgder(1,1,k),auxmat1(1,1)) call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1)) @@ -8482,12 +10581,12 @@ C Cartesian derivatives. s2=0.5d0*scalar2(Ub2(1,i),auxvec(1)) if (j.eq.l+1) then call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat), - & b1(1,itj1),auxvec(1)) - s3=-0.5d0*scalar2(b1(1,itj),auxvec(1)) + & b1(1,j+1),auxvec(1)) + s3=-0.5d0*scalar2(b1(1,j),auxvec(1)) else call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat), - & b1(1,itl1),auxvec(1)) - s3=-0.5d0*scalar2(b1(1,itl),auxvec(1)) + & b1(1,l+1),auxvec(1)) + s3=-0.5d0*scalar2(b1(1,l),auxvec(1)) endif call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1), & pizda(1,1)) @@ -8556,11 +10655,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 @@ -8587,12 +10686,12 @@ cd write (2,*) 'eello6_5',eello6_5 #ifdef MOMENT call transpose2(AEA(1,1,1),auxmat(1,1)) call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1)) - ss1=scalar2(Ub2(1,i+2),b1(1,itl)) + ss1=scalar2(Ub2(1,i+2),b1(1,l)) s1 = (auxmat(1,1)+auxmat(2,2))*ss1 #endif - call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1)) + call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1)) call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1)) - s2 = scalar2(b1(1,itk),vtemp1(1)) + s2 = scalar2(b1(1,k),vtemp1(1)) #ifdef MOMENT call transpose2(AEA(1,1,2),atemp(1,1)) call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1)) @@ -8607,7 +10706,7 @@ cd write (2,*) 'eello6_5',eello6_5 call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1)) call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1)) call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1)) - ss13 = scalar2(b1(1,itk),vtemp4(1)) + ss13 = scalar2(b1(1,k),vtemp4(1)) s13 = (gtemp(1,1)+gtemp(2,2))*ss13 #endif c write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13 @@ -8641,12 +10740,12 @@ C Derivatives in gamma(i+3) #ifdef MOMENT call transpose2(AEA(1,1,1),auxmatd(1,1)) call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) - ss1d=scalar2(Ub2der(1,i+2),b1(1,itl)) + ss1d=scalar2(Ub2der(1,i+2),b1(1,l)) s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d #endif - call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1)) + call matvec2(EUgder(1,1,i+2),b1(1,l),vtemp1d(1)) call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1)) - s2d = scalar2(b1(1,itk),vtemp1d(1)) + s2d = scalar2(b1(1,k),vtemp1d(1)) #ifdef MOMENT call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1)) s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1)) @@ -8694,9 +10793,9 @@ C Derivatives in gamma(i+5) call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 #endif - call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1)) + call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1d(1)) call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1)) - s2d = scalar2(b1(1,itk),vtemp1d(1)) + s2d = scalar2(b1(1,k),vtemp1d(1)) #ifdef MOMENT call transpose2(AEA(1,1,2),atempd(1,1)) call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1)) @@ -8706,7 +10805,7 @@ C Derivatives in gamma(i+5) s12d = scalar2(Ub2(1,i+2),vtemp3d(1)) #ifdef MOMENT call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1)) - ss13d = scalar2(b1(1,itk),vtemp4d(1)) + ss13d = scalar2(b1(1,k),vtemp4d(1)) s13d = (gtemp(1,1)+gtemp(2,2))*ss13d #endif c s1d=0.0d0 @@ -8730,10 +10829,10 @@ C Cartesian derivatives call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1)) s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1 #endif - call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1)) + call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1)) call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1), & vtemp1d(1)) - s2d = scalar2(b1(1,itk),vtemp1d(1)) + s2d = scalar2(b1(1,k),vtemp1d(1)) #ifdef MOMENT call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1)) call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1)) @@ -8777,7 +10876,7 @@ c s13d=0.0d0 derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4), & vtemp4d(1)) - ss13d = scalar2(b1(1,itk),vtemp4d(1)) + ss13d = scalar2(b1(1,k),vtemp4d(1)) s13d = (gtemp(1,1)+gtemp(2,2))*ss13d derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d enddo @@ -9013,4 +11112,701 @@ crc print *,((prod(i,j),i=1,2),j=1,2) 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.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 + return + end +C----------------------------------------------------------------------- +C----------------------------------------------------------- +C This subroutine is to mimic the histone like structure but as well can be +C utilizet to nanostructures (infinit) small modification has to be used to +C make it finite (z gradient at the ends has to be changes as well as the x,y +C gradient has to be modified at the ends +C The energy function is Kihara potential +C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +C 4eps is depth of well sigma is r_minimum r is distance from center of tube +C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +C simple Kihara potential + subroutine calctube(Etube) + 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' + double precision tub_r,vectube(3),enetube(maxres*2) + Etube=0.0d0 + do i=1,2*nres + enetube(i)=0.0d0 + enddo +C first we calculate the distance from tube center +C first sugare-phosphate group for NARES this would be peptide group +C for UNRES + do i=1,nres +C lets ommit dummy atoms for now + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +C now calculate distance from center of tube and direction vectors + vectube(1)=(c(1,i)+c(1,i+1))/2.0d0-tubecenter(1) + vectube(2)=(c(2,i)+c(2,i+1))/2.0d0-tubecenter(2) +C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) + +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6 +C write(iout,*) "TU13",i,rdiff6,enetube(i) +C print *,rdiff,rdiff6,pep_aa_tube +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6+ + & 6.0d0*pep_bb_tube)/rdiff6/rdiff +C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +C &rdiff,fac + +C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + enddo +C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) + do i=1,nres +C Lets not jump over memory as we use many times iti + iti=itype(i) +C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) +C in UNRES uncomment the line below as GLY has no side-chain... +C .or.(iti.eq.10) + & ) cycle + vectube(1)=c(1,i+nres)-tubecenter(1) + vectube(2)=c(2,i+nres)-tubecenter(2) + +C as the tube is infinity we do not calculate the Z-vector use of Z +C as chosen axis + vectube(3)=0.0d0 +C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +C and its 6 power + rdiff6=rdiff**6.0d0 +C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6 +C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +C now we calculate gradient + fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+ + & 6.0d0*sc_bb_tube/rdiff6/rdiff +C now direction of gg_tube vector + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + enddo + do i=1,2*nres + Etube=Etube+enetube(i) + enddo +C print *,"ETUBE", etube + return + end +C TO DO 1) add to total energy +C 2) add to gradient summation +C 3) add reading parameters (AND of course oppening of PARAM file) +C 4) add reading the center of tube +C 5) add COMMONs +C 6) add to zerograd