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=e3a17dd1d69b8645e431f6d8003e0c3c5bda09f5;hb=46f3f544d5a60812e2e11e03b619903f2ed052eb;hp=e9d67bc2c482af499ec9b25402d508c3fe0d4e55;hpb=05c18ccfcefc037cb8394bd779a4e51b7cbea6ec;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 e9d67bc..e3a17dd 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -10,6 +10,8 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include "mpif.h" double precision weights_(n_ene) + integer IERR + integer status(MPI_STATUS_SIZE) #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' @@ -25,6 +27,13 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.CONTROL' include 'COMMON.TIME1' include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' + double precision fac_shieldbuf(maxres), + & grad_shield_locbuf(3,maxcontsshi,-1:maxres), + & grad_shield_sidebuf(3,maxcontsshi,-1:maxres), + & grad_shieldbuf(3,-1:maxres) + integer ishield_listbuf(maxres), + &shield_listbuf(maxcontsshi,maxres) #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -55,6 +64,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) @@ -81,6 +92,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 @@ -146,6 +158,79 @@ C write (iout,*) "shield_mode",shield_mode call set_shield_fac else if (shield_mode.eq.2) then call set_shield_fac2 + if (nfgtasks.gt.1) then +C#define DEBUG +#ifdef DEBUG + write(iout,*) "befor reduce fac_shield reduce" + do i=1,nres + write(2,*) "fac",itype(i),fac_shield(i),grad_shield(1,i) + write(2,*) "list", shield_list(1,i),ishield_list(i), + & grad_shield_side(1,1,i),grad_shield_loc(1,1,i) + enddo +#endif + call MPI_Allgatherv(fac_shield(ivec_start), + & ivec_count(fg_rank1), + & MPI_DOUBLE_PRECISION,fac_shieldbuf(1),ivec_count(0), + & ivec_displ(0), + & MPI_DOUBLE_PRECISION,FG_COMM,IERR) + call MPI_Allgatherv(shield_list(1,ivec_start), + & ivec_count(fg_rank1), + & MPI_I50,shield_listbuf(1,1),ivec_count(0), + & ivec_displ(0), + & MPI_I50,FG_COMM,IERR) + call MPI_Allgatherv(ishield_list(ivec_start), + & ivec_count(fg_rank1), + & MPI_INTEGER,ishield_listbuf(1),ivec_count(0), + & ivec_displ(0), + & MPI_INTEGER,FG_COMM,IERR) + call MPI_Allgatherv(grad_shield(1,ivec_start), + & ivec_count(fg_rank1), + & MPI_UYZ,grad_shieldbuf(1,1),ivec_count(0), + & ivec_displ(0), + & MPI_UYZ,FG_COMM,IERR) + call MPI_Allgatherv(grad_shield_side(1,1,ivec_start), + & ivec_count(fg_rank1), + & MPI_SHI,grad_shield_sidebuf(1,1,1),ivec_count(0), + & ivec_displ(0), + & MPI_SHI,FG_COMM,IERR) + call MPI_Allgatherv(grad_shield_loc(1,1,ivec_start), + & ivec_count(fg_rank1), + & MPI_SHI,grad_shield_locbuf(1,1,1),ivec_count(0), + & ivec_displ(0), + & MPI_SHI,FG_COMM,IERR) + do i=1,nres + fac_shield(i)=fac_shieldbuf(i) + ishield_list(i)=ishield_listbuf(i) + do j=1,3 + grad_shield(j,i)=grad_shieldbuf(j,i) + enddo !j + do j=1,ishield_list(i) + shield_list(j,i)=shield_listbuf(j,i) + do k=1,3 + grad_shield_loc(k,j,i)=grad_shield_locbuf(k,j,i) + grad_shield_side(k,j,i)=grad_shield_sidebuf(k,j,i) + enddo !k + enddo !j + enddo !i +#ifdef DEBUG + write(iout,*) "after reduce fac_shield reduce" + do i=1,nres + write(2,*) "fac",itype(i),fac_shield(i),grad_shield(1,i) + write(2,*) "list", shield_list(1,i),ishield_list(i), + & grad_shield_side(1,1,i),grad_shield_loc(1,1,i) + enddo +#endif +C#undef DEBUG + endif +#ifdef DEBUG + do i=1,nres + write(iout,*) fac_shield(i),ishield_list(i),i,grad_shield(1,i) + do j=1,ishield_list(i) + write(iout,*) "grad", grad_shield_side(1,j,i), + & grad_shield_loc(1,j,i) + enddo + enddo +#endif endif c print *,"Processor",myrank," left VEC_AND_DERIV" if (ipot.lt.6) then @@ -297,6 +382,8 @@ C based on partition function C print *,"przed lipidami" if (wliptran.gt.0) then call Eliptransfer(eliptran) + else + eliptran=0.0d0 endif C print *,"za lipidami" if (AFMlog.gt.0) then @@ -304,6 +391,17 @@ C print *,"za lipidami" else if (selfguide.gt.0) then call AFMvel(Eafmforce) endif + if (TUBElog.eq.1) then +C print *,"just before call" + call calctube(Etube) + elseif (TUBElog.eq.2) then + call calctube2(Etube) + elseif (TUBElog.eq.3) then + call calcnano(Etube) + else + Etube=0.0d0 + endif + #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -348,6 +446,7 @@ C 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" @@ -442,6 +541,7 @@ cMS$ATTRIBUTES C :: proc_proc 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 @@ -449,7 +549,7 @@ cMS$ATTRIBUTES C :: proc_proc & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce - & +ethetacnstr + & +ethetacnstr+wtube*Etube #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc @@ -458,7 +558,7 @@ cMS$ATTRIBUTES C :: proc_proc & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran & +Eafmforce - & +ethetacnstr + & +ethetacnstr+wtube*Etube #endif energia(0)=etot c detecting NaNQ @@ -568,10 +668,31 @@ c enddo & +wturn3*gshieldc_t3(j,i) & +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) + & +wtube*gg_tube(j,i) + enddo - enddo + enddo +C j=1 +C i=0 +C print *,"KUPA2",gradbufc(j,i),wsc*gvdwc(j,i), +C & wscp*gvdwc_scp(j,i),gvdwc_scpp(j,i), +C & welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i), +C & wel_loc*gel_loc_long(j,i), +C & wcorr*gradcorr_long(j,i), +C & wcorr5*gradcorr5_long(j,i), +C & wcorr6*gradcorr6_long(j,i), +C & wturn6*gcorr6_turn_long(j,i), +C & wstrain*ghpbc(j,i) +C & ,wliptran*gliptranc(j,i) +C & ,gradafm(j,i) +C & ,welec*gshieldc(j,i) +C & ,wcorr*gshieldc_ec(j,i) +C & ,wturn3*gshieldc_t3(j,i) +C & ,wturn4*gshieldc_t4(j,i) +C & ,wel_loc*gshieldc_ll(j,i) +C & ,wtube*gg_tube(j,i) #else do i=0,nct do j=1,3 @@ -591,6 +712,8 @@ c enddo & +wcorr*gshieldc_ec(j,i) & +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) + & +wtube*gg_tube(j,i) + enddo @@ -624,7 +747,7 @@ c call flush(iout) #ifdef TIMING c time_allreduce=time_allreduce+MPI_Wtime()-time00 #endif - do i=nnt,nres + do i=0,nres do k=1,3 gradbufc(k,i)=0.0d0 enddo @@ -747,11 +870,7 @@ C print *,gradafm(1,13),"AFM" & +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)+ @@ -784,9 +903,7 @@ C print *,gradafm(1,13),"AFM" & +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 @@ -801,11 +918,47 @@ C print *,gradafm(1,13),"AFM" & +wturn3*gshieldx_t3(j,i) & +wturn4*gshieldx_t4(j,i) & +wel_loc*gshieldx_ll(j,i) + & +wtube*gg_tube_sc(j,i) enddo - enddo + enddo +C i=0 +C j=1 +C print *,"KUPA", gradbufc(j,i),welec*gelc(j,i), +C & wel_loc*gel_loc(j,i), +C & 0.5d0*wscp*gvdwc_scpp(j,i), +C & welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i), +C & wel_loc*gel_loc_long(j,i), +C & wcorr*gradcorr_long(j,i), +C & wcorr5*gradcorr5_long(j,i), +C & wcorr6*gradcorr6_long(j,i), +C & wturn6*gcorr6_turn_long(j,i), +C & wbond*gradb(j,i), +C & wcorr*gradcorr(j,i), +C & wturn3*gcorr3_turn(j,i), +C & wturn4*gcorr4_turn(j,i), +C & wcorr5*gradcorr5(j,i), +C & wcorr6*gradcorr6(j,i), +C & wturn6*gcorr6_turn(j,i), +C & wsccor*gsccorc(j,i) +C & ,wscloc*gscloc(j,i) +C & ,wliptran*gliptranc(j,i) +C & ,gradafm(j,i) +C & +welec*gshieldc(j,i) +C & +welec*gshieldc_loc(j,i) +C & +wcorr*gshieldc_ec(j,i) +C & +wcorr*gshieldc_loc_ec(j,i) +C & +wturn3*gshieldc_t3(j,i) +C & +wturn3*gshieldc_loc_t3(j,i) +C & +wturn4*gshieldc_t4(j,i) +C & ,wturn4*gshieldc_loc_t4(j,i) +C & ,wel_loc*gshieldc_ll(j,i) +C & ,wel_loc*gshieldc_loc_ll(j,i) +C & ,wtube*gg_tube(j,i) + +C print *,gg_tube(1,0),"TU3" #ifdef DEBUG write (iout,*) "gloc before adding corr" do i=1,4*nres @@ -857,7 +1010,7 @@ c#undef DEBUG call MPI_Barrier(FG_COMM,IERR) time_barrier_g=time_barrier_g+MPI_Wtime()-time00 time00=MPI_Wtime() - call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres, + call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) @@ -1101,6 +1254,7 @@ C------------------------------------------------------------------------ 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, @@ -1109,6 +1263,7 @@ C------------------------------------------------------------------------ & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & 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)'/ @@ -1136,6 +1291,7 @@ C------------------------------------------------------------------------ & '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 @@ -1146,6 +1302,7 @@ C------------------------------------------------------------------------ & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, & ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, + & etube,wtube, & etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ @@ -1172,6 +1329,7 @@ C------------------------------------------------------------------------ & '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 @@ -1780,6 +1938,7 @@ C lipbufthick is thickenes of lipid buffore & +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??" @@ -2056,6 +2215,7 @@ C lipbufthick is thickenes of lipid buffore & +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 @@ -2783,28 +2943,28 @@ c write(iout,*) 'nphi=',nphi,nres #endif #ifdef NEWCORR if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp+1 + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp+1 + iti1=nloctyp endif c write(iout,*),i - b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0) + b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0) & +bnew1(2,1,iti)*dsin(theta(i-1)) - & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0) + & +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0) gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 & +bnew1(2,1,iti)*dcos(theta(i-1)) & -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0 c & +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) - b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0) + b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0) & +bnew2(2,1,iti)*dsin(theta(i-1)) - & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0) + & +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0) c & +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i)) c &*(cos(theta(i)/2.0) gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0 @@ -2843,15 +3003,15 @@ c write (iout,*) 'theta=', theta(i-1) enddo #else 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 b1(1,i-2)=b(3,iti) b1(2,i-2)=b(5,iti) @@ -2940,15 +3100,15 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then endif c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itype2loc(itype(i-2)) else - iti=ntortyp + iti=nloctyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp + iti1=nloctyp endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) @@ -2961,8 +3121,8 @@ c if (i .gt. iatel_s+2) then call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" #endif -c write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti), -c & EE(1,2,iti),EE(2,2,iti) +c write(iout,*) "co jest kurwa", iti, EE(1,1,i),EE(2,1,i), +c & EE(1,2,iti),EE(2,2,i) call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) c write(iout,*) "Macierz EUG", @@ -2997,18 +3157,24 @@ c & eug(2,2,i-2) c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then if (itype(i-1).le.ntyp) then - iti1 = itortyp(itype(i-1)) + iti1 = itype2loc(itype(i-1)) else - iti1=ntortyp + iti1=nloctyp endif else - iti1=ntortyp + iti1=nloctyp endif do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -C write (iout,*) 'mumu',i,b1(1,i-1),Ub2(1,i-2) -c write (iout,*) 'mu ',mu(:,i-2),i-2 +#ifdef MUOUT + write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1), + & rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2), + & -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2), + & dsqrt(b2(1,i-1)**2+b2(2,i-1)**2) + & +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2), + & ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itype2loc(iti)) +#endif cd write (iout,*) '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) @@ -3295,11 +3461,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 @@ -3333,6 +3499,7 @@ C 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), @@ -3417,21 +3584,23 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end - if (i.le.1) cycle +c if (i.le.1) cycle C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+4).gt.nres) - & .or.((i-1).le.0) +C Adam: Unnecessary: handled by iturn3_end and iturn3_start +c & .or.((i+4).gt.nres) +c & .or.((i-1).le.0) C end of changes by Ana & .or. itype(i+2).eq.ntyp1 & .or. itype(i+3).eq.ntyp1) cycle - if(i.gt.1)then - if(itype(i-1).eq.ntyp1)cycle - end if - if(i.LT.nres-3)then - if (itype(i+4).eq.ntyp1) cycle - end if +C Adam: Instructions below will switch off existing interactions +c if(i.gt.1)then +c if(itype(i-1).eq.ntyp1)cycle +c end if +c if(i.LT.nres-3)then +c if (itype(i+4).eq.ntyp1) cycle +c end if dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3447,23 +3616,47 @@ C end of changes by Ana if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + zmedi2=mod(zmedi,boxzsize) + if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize + if ((zmedi2.gt.bordlipbot) + &.and.(zmedi2.lt.bordliptop)) then +C the energy transfer exist + if (zmedi2.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zmedi2-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi2.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0d0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0d0 + endif 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.le.1) cycle + 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 - & .or.((i+5).gt.nres) - & .or.((i-1).le.0) +c & .or.((i+5).gt.nres) +c & .or.((i-1).le.0) C end of changes suggested by Ana & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +c & .or. itype(i+5).eq.ntyp1 +c & .or. itype(i).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3505,7 +3698,30 @@ c endif if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize - + zmedi2=mod(zmedi,boxzsize) + if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize + if ((zmedi2.gt.bordlipbot) + &.and.(zmedi2.lt.bordliptop)) then +C the energy transfer exist + if (zmedi2.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zmedi2-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi2.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif num_conti=num_cont_hb(i) c write(iout,*) "JESTEM W PETLI" call eelecij(i,i+3,ees,evdw1,eel_loc) @@ -3523,14 +3739,14 @@ c CTU KURWA do i=iatel_s,iatel_e C do i=75,75 - if (i.le.1) cycle +c if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds - & .or.((i+2).gt.nres) - & .or.((i-1).le.0) +c & .or.((i+2).gt.nres) +c & .or.((i-1).le.0) C end of changes by Ana - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +c & .or. itype(i+2).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -3547,6 +3763,29 @@ C end of changes by Ana if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + if ((zmedi.gt.bordlipbot) + &.and.(zmedi.lt.bordliptop)) then +C the energy transfer exist + if (zmedi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zmedi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi)/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 print *,sslipi,"TU?!" C xmedi=xmedi+xshift*boxxsize C ymedi=ymedi+yshift*boxysize C zmedi=zmedi+zshift*boxzsize @@ -3583,14 +3822,14 @@ C I TU KURWA do j=ielstart(i),ielend(i) C do j=16,17 C write (iout,*) i,j - if (j.le.1) cycle +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 - & .or.((j+2).gt.nres) - & .or.((j-1).le.0) +c & .or.((j+2).gt.nres) +c & .or.((j-1).le.0) C end of changes by Ana - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +c & .or.itype(j+2).eq.ntyp1 +c & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j @@ -3651,6 +3890,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 @@ -3680,6 +3920,28 @@ C zj=c(3,j)+0.5D0*dzj-zmedi 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" + 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 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj @@ -3774,13 +4036,20 @@ C fac_shield(j)=0.6 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) ees=ees+eesij +C FOR NOW SHIELD IS NOT USED WITH LIPSCALE +C & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) ees=ees+eesij + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) +C print *,"TUCC",(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 endif evdw1=evdw1+evdwij*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) +C print *,sslipi,sslipj,lipscale**2, +C & (sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 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, @@ -3790,6 +4059,7 @@ cd & xmedi,ymedi,zmedi,xj,yj,zj write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &'evdw1',i,j,evdwij &,iteli,itelj,aaa,evdw1 + write (iout,*) sss write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, &fac_shield(i),fac_shield(j) endif @@ -3799,7 +4069,9 @@ C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE facvdw=-6*rrmij*(ev1+evdwij)*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) facel=-3*rrmij*(el1+eesij) + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) fac1=fac erij(1)=xj*rmij erij(2)=yj*rmij @@ -3898,6 +4170,12 @@ 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) +C Lipidic part for lipscale + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2 + + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. @@ -3909,8 +4187,13 @@ cgrad enddo cgrad enddo if (sss.gt.0.0) then ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) else ggg(1)=0.0 ggg(2)=0.0 @@ -3926,6 +4209,12 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 + * * Loop over residues i+1 thru j-1. * @@ -3937,6 +4226,7 @@ cgrad enddo #else C MARYSIA facvdw=(ev1+evdwij)*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -3972,12 +4262,22 @@ cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 + #endif * * Angular part @@ -3996,6 +4296,7 @@ cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* & fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -4017,10 +4318,12 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) & +((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 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) 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)) & *fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -4235,6 +4538,8 @@ C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + C Now derivative over eel_loc if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then @@ -4291,6 +4596,8 @@ C Calculate patrial derivative for theta angle & +a32*gmuij1(3) & +a33*gmuij1(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4307,6 +4614,8 @@ c & a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4320,6 +4629,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) geel_loc_ji= & +a22*gmuji2(1) @@ -4332,11 +4642,13 @@ 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + #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 + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2f7.3)') + & 'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2) 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) @@ -4348,22 +4660,35 @@ C Partial derivatives in virtual-bond dihedral angles gamma & (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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + 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) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/ + & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/ + & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -4374,18 +4699,22 @@ C Remaining derivatives of eello 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo ENDIF @@ -4631,7 +4960,42 @@ C Third- and fourth-order contributions from turns & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 j=i+2 -c write (iout,*) "eturn3",i,j,j1,j2 +C xj=(c(1,j)+c(1,j+1))/2.0d0 +C yj=(c(2,j)+c(2,j+1))/2.0d0 + zj=(c(3,j)+c(3,j+1))/2.0d0 +C xj=mod(xj,boxxsize) +C if (xj.lt.0) xj=xj+boxxsize +C yj=mod(yj,boxysize) +C 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" + 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 +C sslipj=0.0 +C ssgradlipj=0.0d0 + +C write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -4660,24 +5024,35 @@ c auxalary matrix for i+2 and constant i+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 + fac_shield(i)=1.0d0 + fac_shield(j)=1.0d0 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)) +C if (j.eq.78) +C & write(iout,*) i,j,fac_shield(i),fac_shield(j) + eello_turn3=eello_turn3+ +C & 1.0*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + &0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) - eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t3= + &0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) +#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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) +#endif C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') C Derivatives in shield mode @@ -4733,6 +5108,8 @@ C Derivatives in gamma(i) 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + 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)) @@ -4740,6 +5117,8 @@ C Derivatives in gamma(i+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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -4754,6 +5133,7 @@ c ghalf4=0.5d0*agg(l,4) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -4763,6 +5143,7 @@ c ghalf4=0.5d0*agg(l,4) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -4771,6 +5152,8 @@ c ghalf4=0.5d0*agg(l,4) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -4779,7 +5162,18 @@ c ghalf4=0.5d0*agg(l,4) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo + gshieldc_t3(3,i)=gshieldc_t3(3,i)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j)=gshieldc_t3(3,j)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + +C print *,ssgradlipi,ssgradlipj,eello_t3,sslipi,sslipj return end C------------------------------------------------------------------------------- @@ -4828,13 +5222,44 @@ 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" + zj=(c(3,j)+c(3,j+1))/2.0d0 +C xj=mod(xj,boxxsize) +C if (xj.lt.0) xj=xj+boxxsize +C yj=mod(yj,boxysize) +C if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize +C if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" + 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 + 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)) @@ -4909,6 +5334,8 @@ C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -4968,13 +5395,17 @@ cd & ' eello_turn4_num',8*eello_turn4_num gloc(nphi+i,icg)=gloc(nphi+i,icg) & -(gs13+gsE13+gsEE1)*wturn4 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg) & -(gs23+gs21+gsEE2)*wturn4 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg) & -(gs32+gsE31+gsEE3)*wturn4 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- c & gs2 @@ -4992,6 +5423,8 @@ C Derivatives in gamma(i) 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + 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)) @@ -5001,6 +5434,8 @@ C Derivatives in gamma(i+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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + 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)) @@ -5013,6 +5448,8 @@ C Derivatives in gamma(i+2) 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -5033,6 +5470,8 @@ C Derivatives of this turn contributions in DC(i+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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + enddo endif C Remaining derivatives of this turn contribution @@ -5052,6 +5491,8 @@ C Remaining derivatives of this turn contribution 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -5067,6 +5508,8 @@ C Remaining derivatives of this turn contribution 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -5082,6 +5525,8 @@ C Remaining derivatives of this turn contribution 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -5098,7 +5543,16 @@ C Remaining derivatives of this turn contribution 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) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo + gshieldc_t4(3,i)=gshieldc_t4(3,i)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j)=gshieldc_t4(3,j)+ + & ssgradlipj*eello_t4/4.0d0*lipscale + gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ + & ssgradlipj*eello_t4/4.0d0*lipscale return end C----------------------------------------------------------------------------- @@ -5840,6 +6294,7 @@ 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 @@ -5853,6 +6308,7 @@ C NO vbldp0 is the equlibrium lenght of spring for peptide group c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) 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 @@ -5873,6 +6329,9 @@ c else do j=1,nbi diff=vbld(i+nres)-vbldsc0(j,iti) + if (energy_dec) write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(j,iti),diff, + & AKSC(j,iti),AKSC(j,iti)*diff*diff ud(j)=aksc(j,iti)*diff u(j)=abond0(j,iti)+0.5d0*ud(j)*diff enddo @@ -6206,6 +6665,7 @@ C etheta=0.0D0 do i=ithet_start,ithet_end c print *,i,itype(i-1),itype(i),itype(i-2) +C 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 print *,i,theta(i) @@ -7473,11 +7933,12 @@ C The rigorous attempt to derive energy function include 'COMMON.TORCNSTR' include 'COMMON.CONTROL' logical lprn - double precision thybt1(maxtermkcc),thybt2(maxtermkcc) +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 @@ -7497,60 +7958,86 @@ c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle sumnonchebyshev=0.0d0 sumchebyshev=0.0d0 C to avoid multiple devision by 2 - theti22=0.5d0*theta(i) +c theti22=0.5d0*theta(i) C theta 12 is the theta_1 /2 C theta 22 is theta_2 /2 - theti12=0.5d0*theta(i-1) +c theti12=0.5d0*theta(i-1) C and appropriate sinus function - sinthet2=dsin(theta(i)) 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 + 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) - sint1t2n=sint1t2**j - sumnonchebyshev= - & sint1t2n*(v1ij*cosphi+v2ij*sinphi) - actval=sint1t2n*(v1ij*cosphi+v2ij*sinphi) + 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 - glocig=j*(v2ij*cosphi-v1ij*sinphi)*sint1t2n + 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=actval/sinthet1*j*costhet1 - glocit2=actval/sinthet2*j*costhet2 + 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 - do k=1,nterm_kcc_Tb(itori,itori1) - thybt1(k)=v1_chyb(k,j,itori,itori1) - thybt2(k)=v2_chyb(k,j,itori,itori1) +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 - enddo - sumth1thyb=tschebyshev - & (1,nterm_kcc_Tb(itori,itori1),thybt1(1),dcos(theti12)**2) - gradthybt1=gradtschebyshev - & (0,nterm_kcc_Tb(itori,itori1)-1,thybt1(1), - & dcos(theti12)**2) - & *dcos(theti12)*(-dsin(theti12)) - sumth2thyb=tschebyshev - & (1,nterm_kcc_Tb(itori,itori1),thybt2(1),dcos(theti22)**2) - gradthybt2=gradtschebyshev - & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1), - & dcos(theti22)**2) - & *dcos(theti22)*(-dsin(theti22)) +c enddo C print *, sumth1thyb, gradthybt1, sumth2thyb, gradthybt2, C & gradtschebyshev C & (0,nterm_kcc_Tb(itori,itori1)-1,thybt2(1), @@ -7558,22 +8045,19 @@ C & dcos(theti22)**2), C & dsin(theti22) C now overal sumation - etors=etors+(1.0d0+sumth1thyb+sumth2thyb)*sumnonchebyshev 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 - & *(1.0d0+sumth1thyb+sumth2thyb) + 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*(1.0d0+sumth1thyb+sumth2thyb)+ - & sumnonchebyshev*gradthybt1) + 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*(1.0d0+sumth1thyb+sumth2thyb)+ - & sumnonchebyshev*gradthybt2) - enddo + 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 @@ -7622,7 +8106,8 @@ C Set lprn=.true. for debugging lprn=.false. c lprn=.true. C print *,"wchodze kcc" - if (tormode.ne.2) etheta=0.0D0 + 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 @@ -7635,6 +8120,8 @@ c print *,i,itype(i-1),itype(i),itype(i-2) 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) @@ -7643,7 +8130,7 @@ C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0) gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang* & gradthybt1*sinthet*(-0.5d0) enddo - if (tormode.ne.2) then + if (tor_mode.ne.2) then ethetacnstr=0.0d0 C print *,ithetaconstr_start,ithetaconstr_end,"TU" do i=ithetaconstr_start,ithetaconstr_end @@ -7736,6 +8223,9 @@ c 3 = SC...Ca...Ca...SCi esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo + if (energy_dec) write(iout,'(a9,2i4,f8.3,3i4)') "esccor",i,j, + & esccor,intertyp, + & isccori, isccori1 c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci if (lprn) @@ -8839,9 +9329,9 @@ C--------------------------------------------------------------------------- & auxmat(2,2) iti1 = itortyp(itype(i+1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1)) + itj1 = itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -8929,16 +9419,16 @@ cd write (iout,*) "a_chujkl",((a_chuj(iii,jjj,kk,k),iii=1,2),jjj=1,2) if (l.eq.j+1) then C parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itype2loc(itype(i)) else - iti=ntortyp + iti=nloctyp endif - itk1=itortyp(itype(k+1)) - itj=itortyp(itype(j)) + itk1=itype2loc(itype(k+1)) + itj=itype2loc(itype(j)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -9082,17 +9572,17 @@ C End vectors else C Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itype2loc(itype(i)) else - iti=ntortyp + iti=nloctyp endif - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + itk1=itype2loc(itype(k+1)) + itl=itype2loc(itype(l)) + itj=itype2loc(itype(j)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif C A2 kernel(j-1)T A1T call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), @@ -9430,9 +9920,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 @@ -9501,7 +9991,7 @@ C Cartesian gradient c goto 1112 c1111 continue C Contribution from graph II - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9582,7 +10072,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV cd1110 continue - call transpose2(EE(1,1,itl),auxmat(1,1)) + call transpose2(EE(1,1,l),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9655,7 +10145,7 @@ C Cartesian gradient cd goto 1112 C Contribution from graph IV 1110 continue - call transpose2(EE(1,1,itj),auxmat(1,1)) + call transpose2(EE(1,1,j),auxmat(1,1)) call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -9952,7 +10442,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)) @@ -10244,16 +10734,16 @@ C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. iti=itortyp(itype(i)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) + itk=itype2loc(itype(k)) + itk1=itype2loc(itype(k+1)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) @@ -10262,7 +10752,7 @@ C energy moment and not to the cluster cumulant. s2=0.5d0*scalar2(b1(1,k),auxvec(1)) call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1)) s3=0.5d0*scalar2(b1(1,j+1),auxvec(1)) - call transpose2(EE(1,1,itk),auxmat(1,1)) + call transpose2(EE(1,1,k),auxmat(1,1)) call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1)) vv(1)=pizda(1,1)+pizda(2,2) vv(2)=pizda(2,1)-pizda(1,2) @@ -10360,24 +10850,24 @@ C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. cd write (2,*) 'eello_graph4: wturn6',wturn6 - iti=itortyp(itype(i)) - itj=itortyp(itype(j)) + iti=itype2loc(itype(i)) + itj=itype2loc(itype(j)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itype2loc(itype(j+1)) else - itj1=ntortyp + itj1=nloctyp endif - itk=itortyp(itype(k)) + itk=itype2loc(itype(k)) if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) + itk1=itype2loc(itype(k+1)) else - itk1=ntortyp + itk1=nloctyp endif - itl=itortyp(itype(l)) + itl=itype2loc(itype(l)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itype2loc(itype(l+1)) else - itl1=ntortyp + itl1=nloctyp endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, @@ -10597,11 +11087,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 @@ -11085,7 +11575,7 @@ C do i=1,1 if (itype(i).eq.ntyp1) cycle positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) - if (positi.le.0) positi=positi+boxzsize + 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 @@ -11441,7 +11931,7 @@ C-------------------------------------------------------------------------- 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 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 @@ -11467,18 +11957,29 @@ C first for shielding is setting of function of side-chains include 'COMMON.IOUNITS' include 'COMMON.SHIELD' include 'COMMON.INTERACT' + include 'COMMON.LOCAL' + 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 write(2,*) "ivec",ivec_start,ivec_end + do i=1,nres + fac_shield(i)=0.0d0 + do j=1,3 + grad_shield(j,i)=0.0d0 + enddo + enddo 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 + do i=ivec_start,ivec_end +C do i=1,nres-1 +C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle ishield_list(i)=0 + if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle 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 @@ -11511,6 +12012,7 @@ C now sscale fraction C print *,buff_shield,"buff" C now sscale if (sh_frac_dist.le.0.0) cycle +C print *,ishield_list(i),i 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 @@ -11625,8 +12127,696 @@ C grad_shield_side is Cbeta sidechain gradient 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) +C write(2,*) "TOTAL VOLUME",i,itype(i),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=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=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=itube_start,itube_end +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 + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-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) +C print *,gg_tube(1,0),"TU" + + + do i=itube_start,itube_end +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 + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp +C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +C & tubecenter(2) + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-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=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + 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 + +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 calctube2(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=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=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=itube_start,itube_end +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 +C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) +C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize +C vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) +C if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-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 THIS FRAGMENT MAKES TUBE FINITE + positi=mod((c(3,i)+c(3,i+1))/2.0d0,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) + print *,positi,bordtubebot,buftubebot,bordtubetop + if ((positi.gt.bordtubebot) + & .and.(positi.lt.bordtubetop)) then +C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- + & ((positi-bordtubebot)/tubebufthick) +C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick + print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i)=enetube(i)+sstube*tubetranenepep +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- + &((bordtubetop-positi)/tubebufthick) + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i)=enetube(i)+sstube*tubetranenepep +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i)=enetube(i)+sstube*tubetranenepep +C print *,"I am in true lipid" + endif + else +C sstube=0.0d0 +C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor + +C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=enetube(i)+sstube* + &(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*sstube +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 + gg_tube(3,i)=gg_tube(3,i) + &+ssgradtube*enetube(i)/sstube/2.0d0 + gg_tube(3,i-1)= gg_tube(3,i-1) + &+ssgradtube*enetube(i)/sstube/2.0d0 + + enddo +C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) +C print *,gg_tube(1,0),"TU" + do i=itube_start,itube_end +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... + & .or.(iti.eq.10) + & ) cycle + vectube(1)=c(1,i+nres) + vectube(1)=mod(vectube(1),boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=c(2,i+nres) + vectube(2)=mod(vectube(2),boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) +C THIS FRAGMENT MAKES TUBE FINITE + 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) + print *,positi,bordtubebot,buftubebot,bordtubetop + if ((positi.gt.bordtubebot) + & .and.(positi.lt.bordtubetop)) then +C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- + & ((positi-bordtubebot)/tubebufthick) +C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick + print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- + &((bordtubetop-positi)/tubebufthick) + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +C print *,"I am in true lipid" + endif + else +C sstube=0.0d0 +C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor +CEND OF FINITE FRAGMENT +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) + & *sstube+enetube(i+nres) +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)*sstube +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 + gg_tube_SC(3,i)=gg_tube_SC(3,i) + &+ssgradtube*enetube(i+nres)/sstube + gg_tube(3,i-1)= gg_tube(3,i-1) + &+ssgradtube*enetube(i+nres)/sstube + + enddo + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + 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 + + +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 calcnano(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), + & enecavtube(maxres*2) + Etube=0.0d0 + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=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=itube_start,itube_end +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 + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + + do j=-1,1 + vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=dabs(vectube(1)-tubecenter(1)) + yminact=dabs(vectube(2)-tubecenter(2)) + zminact=dabs(vectube(3)-tubecenter(3)) + + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + if (zmin.gt.zminact) then + zmin=zminact + ztemp=vectube(3) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(3)=ztemp + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + vectube(3)=vectube(3)-tubecenter(3) + +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 +C 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 + vectube(3)=vectube(3)/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 + if (acavtubpep.eq.0.0d0) then +C go to 667 + enecavtube(i)=0.0 + faccav=0.0 + else + denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6) + enecavtube(i)= + & (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep) + & /denominator + enecavtube(i)=0.0 + faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff)) + & *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff) + & +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) + & /denominator**2.0d0 +C faccav=0.0 +C fac=fac+faccav +C 667 continue + endif +C print *,"TUT",i,iti,rdiff,rdiff6,acavtubpep,denominator, +C & enecavtube(i),faccav +C print *,"licz=", +C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) +CX print *,"finene=",enetube(i+nres)+enecavtube(i) + +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 + + do i=itube_start,itube_end + enecavtube(i)=0.0d0 +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 + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + do j=-1,1 + vectube(1)=dmod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=dmod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + vectube(3)=dmod((c(3,i+nres)),boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=dabs(vectube(1)-tubecenter(1)) + yminact=dabs(vectube(2)-tubecenter(2)) + zminact=dabs(vectube(3)-tubecenter(3)) + + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + if (zmin.gt.zminact) then + zmin=zminact + ztemp=vectube(3) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(3)=ztemp + +C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +C & tubecenter(2) + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + vectube(3)=vectube(3)-tubecenter(3) +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 + vectube(3)=vectube(3)/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 enetube(i+nres)=0.0d0 +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 fac=0.0 +C now direction of gg_tube vector +C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) + if (acavtub(iti).eq.0.0d0) then +C go to 667 + enecavtube(i+nres)=0.0d0 + faccav=0.0d0 + else + denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6) + enecavtube(i+nres)= + & (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti)) + & /denominator +C enecavtube(i)=0.0 + faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff)) + & *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff) + & +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0) + & /denominator**2.0d0 +C faccav=0.0 + fac=fac+faccav +C 667 continue + endif +C print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator, +C & enecavtube(i),faccav +C print *,"licz=", +C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) +C print *,"finene=",enetube(i+nres)+enecavtube(i) + 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 +C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) +C do i=itube_start,itube_end +C enecav(i)=0.0 +C iti=itype(i) +C if (acavtub(iti).eq.0.0) cycle + + + + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) + & +enecavtube(i+nres) + 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