#ifdef MPI
include "mpif.h"
double precision weights_(n_ene)
+ integer IERR
+ integer status(MPI_STATUS_SIZE)
#endif
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
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
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)
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
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.gt.0) then
+ if (shield_mode.eq.1) then
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
C Calculate the virtual-bond-angle energy.
C
if (wang.gt.0d0) then
+ if ((tor_mode.eq.0).or.(tor_mode.eq.2)) then
call ebend(ebe,ethetacnstr)
+ endif
+C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+C energy function
+ if ((tor_mode.eq.1).or.(tor_mode.eq.2)) then
+ call ebend_kcc(ebe,ethetacnstr)
+ endif
else
ebe=0
ethetacnstr=0
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
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
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
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
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"
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
& +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
& +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
& +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
+ 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
& +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
#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
& +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)+
& +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
& +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
+ 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
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)
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
#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
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,wsccro,edihcnstr,
+ & 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)'/
& '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
& 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)'/
& '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
& +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??"
& +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
#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
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)
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)
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",
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)
#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
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),
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)
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)
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)
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)
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
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
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
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
if (shield_mode.gt.0) then
C fac_shield(i)=0.4
C fac_shield(j)=0.6
- el1=el1*fac_shield(i)*fac_shield(j)
- el2=el2*fac_shield(i)*fac_shield(j)
+ 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
+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,
write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
&'evdw1',i,j,evdwij
&,iteli,itelj,aaa,evdw1
- write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+ write (iout,*) sss
+ write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
+ &fac_shield(i),fac_shield(j)
endif
C
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
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)
+ & +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)
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)
+ & +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)
do k=1,3
gshieldc(k,i)=gshieldc(k,i)+
- & grad_shield(k,i)*eesij/fac_shield(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)
+ & 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)
+ & 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)
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
enddo
endif
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.
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
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.
*
#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)
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
cd & (dcosg(k),k=1,3)
do k=1,3
ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
- & fac_shield(i)*fac_shield(j)
+ & 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)
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))
- & *fac_shield(i)*fac_shield(j)
+ & *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)*fac_shield(j)
+ & *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
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)
+ &*((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
+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)
+ geel_loc_ij=(a22*gmuij1(1)
& +a23*gmuij1(2)
& +a32*gmuij1(3)
- & +a33*gmuij1(4)
+ & +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)
& +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)
& +a23*gmuji1(2)
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)
& +a23*gmuji2(2)
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
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)
+ &*((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)
+ & (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)
+ 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)
do l=1,3
gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
& aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
+ & *fac_shield(i)*fac_shield(j)
+ &*((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
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
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
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
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),
& 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
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))
- eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+ if (shield_mode.eq.0) then
+ 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
+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)
+ &*((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)
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
+#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
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)
+ &*((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))
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)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
C Cartesian derivatives
do l=1,3
c ghalf1=0.5d0*agg(l,1)
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)
+ &*((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)
a_temp(2,1)=aggi1(l,3)!+agg(l,3)
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)
+ &*((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
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)
+ &*((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)
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)
+ &*((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-------------------------------------------------------------------------------
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),
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))
-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))
-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))
+ 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=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))
+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
gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
-
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.6
+C fac_shield(j)=0.4
+ endif
eello_turn4=eello_turn4-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+ &*((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)
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)
+ &*((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
#endif
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)
+ &*((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))
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)
+ &*((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))
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)
+ &*((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
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)
+ &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
enddo
endif
C Remaining derivatives of this turn contribution
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+ &*((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)
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)
+ &*((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)
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)
+ &*((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)
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)
+ &*((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-----------------------------------------------------------------------------
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
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
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
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)
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
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)
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.SHIELD'
double precision gx(3),gx1(3)
logical lprn
lprn=.false.
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
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)
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)
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
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
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
& 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)
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
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),
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
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)
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)
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)
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))
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)
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)
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,
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
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
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
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'
+ 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=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
+ 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 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
+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,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)=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
+ vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ zminact=abs(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.0+dcavtubpep*rdiff6*rdiff6)
+ enecavtube(i)=
+ & (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep)
+ & /denominator
+ enecavtube(i)=0.0
+ faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff))
+ & *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(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.0
+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)=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
+ vectube(3)=mod((c(3,i+nres)),boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ zminact=abs(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.0
+ faccav=0.0
+ else
+ denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
+ enecavtube(i+nres)=
+ & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+ & /denominator
+C enecavtube(i)=0.0
+ faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff))
+ & *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(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
projects / unres.git / blobdiff