c if (dyn_ss) call dyn_set_nss
c print *,"Processor",myrank," computed USCSC"
+c write (iout,*) "SCSC computed OK"
+c call flush_(iout)
#ifdef TIMING
time01=MPI_Wtime()
#endif
C the shielding factor is set this factor is describing how each
C peptide group is shielded by side-chains
C the matrix - shield_fac(i) the i index describe the ith between i and i+1
- write (iout,*) "shield_mode",shield_mode
+C write (iout,*) "shield_mode",shield_mode
if (shield_mode.gt.0) then
call set_shield_fac
endif
c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
c & eello_turn4)
endif
+c write (iout,*) "eelec computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed UELEC"
C
C Calculate excluded-volume interaction energy between peptide groups
c write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
+c write (iout,*) "escp computed OK"
+c call flush_(iout)
c
c Calculate the bond-stretching energy
c
call ebond(estr)
+c write (iout,*) "ebond computed OK"
+c call flush_(iout)
C
C Calculate the disulfide-bridge and other energy and the contributions
C from other distance constraints.
ebe=0
ethetacnstr=0
endif
+c write (iout,*) "ebend computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed UB"
C
C Calculate the SC local energy.
C
C print *,"TU DOCHODZE?"
call esc(escloc)
+c write (iout,*) "esc computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed USC"
C
C Calculate the virtual-bond torsional energy.
etors=0
edihcnstr=0
endif
+c write (iout,*) "etor computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed Utor"
C
C 6/23/01 Calculate double-torsional energy
else
etors_d=0
endif
+c write (iout,*) "etor_d computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed Utord"
C
C 21/5/07 Calculate local sicdechain correlation energy
else
esccor=0.0d0
endif
+c write (iout,*) "eback_sc_corr computed OK"
+c call flush_(iout)
C print *,"PRZED MULIt"
c print *,"Processor",myrank," computed Usccorr"
C
n_corr1=0
if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0
& .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then
+c write (iout,*) "Calling multibody_eello"
+c call flush_(iout)
call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
-cd write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
-cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+c write(iout,*)
+c & 'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
+c & " ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+c call flush_(iout)
else
ecorr=0.0d0
ecorr5=0.0d0
eturn6=0.0d0
endif
if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
+c write (iout,*) "Calling multibody_gb_ecorr"
+c call flush_(iout)
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
-cd write (iout,*) "multibody_hb ecorr",ecorr
+c write (iout,*) "Exited multibody_hb ecorr",ecorr
+c call flush_(iout)
endif
+c write (iout,*) "multibody computed OK"
+c call flush_(iout)
c print *,"Processor",myrank," computed Ucorr"
C
C If performing constraint dynamics, call the constraint energy
if (wliptran.gt.0) then
call Eliptransfer(eliptran)
endif
-C print *,"za lipidami"
+c write (iout,*) "lipid energy computed OK"
+c call flush_(iout)
if (AFMlog.gt.0) then
call AFMforce(Eafmforce)
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
endif
+c write (iout,*) "AFMforce computed OK"
+c call flush_(iout)
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
c per molecule then we sum it up in sum_energy subroutine
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
+c write (iout,*) "sum energy OK"
+c call flush_(iout)
if (dyn_ss) call dyn_set_nss
+c write (iout,*) "Exiting energy"
+c call flush_(iout)
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
time_sumene=time_sumene+MPI_Wtime()-time00
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
enddo
enddo
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
enddo
enddo
do i=-1,nct
do j=1,3
#ifdef SPLITELE
+C print *,gradbufc(1,13)
+C print *,welec*gelc(1,13)
+C print *,wel_loc*gel_loc(1,13)
+C print *,0.5d0*(wscp*gvdwc_scpp(1,13))
+C print *,welec*gelc_long(1,13)+wvdwpp*gvdwpp(1,13)
+C print *,wel_loc*gel_loc_long(1,13)
+C print *,gradafm(1,13),"AFM"
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& 0.5d0*(wscp*gvdwc_scpp(j,i)+
& +wscloc*gscloc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
+ & +welec*gshieldc_loc(j,i)
+
+
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
& wel_loc*gel_loc(j,i)+
& +wscloc*gscloc(j,i)
& +wliptran*gliptranc(j,i)
& +gradafm(j,i)
+ & +welec*gshieldc(j,i)
+ & +welec*gshieldc_loc(j,i)
+
#endif
gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
& wsccor*gsccorx(j,i)
& +wscloc*gsclocx(j,i)
& +wliptran*gliptranx(j,i)
+ & +welec*gshieldx(j,i)
enddo
enddo
#ifdef DEBUG
c write (iout,*) 'theta=', theta(i-1)
enddo
#else
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itortyp(itype(i-2))
+ else
+ iti=ntortyp+1
+ 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))
+ else
+ iti1=ntortyp+1
+ endif
b1(1,i-2)=b(3,iti)
b1(2,i-2)=b(5,iti)
b2(1,i-2)=b(2,iti)
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
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
+CTU KURWA
do i=iatel_s,iatel_e
+C do i=75,75
if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
+C I TU KURWA
do j=ielstart(i),ielend(i)
+C do j=16,17
C write (iout,*) i,j
if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
el2=fac4*fac
C MARYSIA
- eesij=(el1+el2)
+C eesij=(el1+el2)
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
if (shield_mode.gt.0) then
- ees=ees+eesij*fac_shield(i)*fac_shield(j)
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ el1=el1*fac_shield(i)*fac_shield(j)
+ el2=el2*fac_shield(i)*fac_shield(j)
+ eesij=(el1+el2)
+ ees=ees+eesij
else
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ eesij=(el1+el2)
ees=ees+eesij
endif
evdw1=evdw1+evdwij*sss
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
+
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
ggg(1)=facel*xj
ggg(2)=facel*yj
ggg(3)=facel*zj
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C if (iresshield.gt.i) then
+C do ishi=i+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C else
+C do ishi=iresshield,i
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C endif
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C if (iresshield.gt.j) then
+C do ishi=j+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C
+C enddo
+C else
+C do ishi=iresshield,j
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
+C enddo
+C endif
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,i)+
+ & grad_shield(k,i)*eesij/fac_shield(i)
+ gshieldc(k,j)=gshieldc(k,j)+
+ & grad_shield(k,j)*eesij/fac_shield(j)
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & grad_shield(k,i)*eesij/fac_shield(i)
+ gshieldc(k,j-1)=gshieldc(k,j-1)+
+ & grad_shield(k,j)*eesij/fac_shield(j)
+
+ enddo
+ endif
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
c gelc(k,j)=gelc(k,j)+ghalf
c enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
+C print *,"before", gelc_long(1,i), gelc_long(1,j)
do k=1,3
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+C & +grad_shield(k,j)*eesij/fac_shield(j)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+C & +grad_shield(k,i)*eesij/fac_shield(i)
+C gelc_long(k,i-1)=gelc_long(k,i-1)
+C & +grad_shield(k,i)*eesij/fac_shield(i)
+C gelc_long(k,j-1)=gelc_long(k,j-1)
+C & +grad_shield(k,j)*eesij/fac_shield(j)
enddo
+C print *,"bafter", gelc_long(1,i), gelc_long(1,j)
+
*
* Loop over residues i+1 thru j-1.
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
ggg(1)=fac*xj
+C+eesij*grad_shield(1,i)+eesij*grad_shield(1,j)
ggg(2)=fac*yj
+C+eesij*grad_shield(2,i)+eesij*grad_shield(2,j)
ggg(3)=fac*zj
+C+eesij*grad_shield(3,i)+eesij*grad_shield(3,j)
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
cd & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
+ & fac_shield(i)*fac_shield(j)
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
+C print *,"before22", gelc_long(1,i), gelc_long(1,j)
do k=1,3
gelc(k,i)=gelc(k,i)
- & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+ & *fac_shield(i)*fac_shield(j)
gelc(k,j)=gelc(k,j)
- & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+ & *fac_shield(i)*fac_shield(j)
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+C print *,"before33", gelc_long(1,i), gelc_long(1,j)
+
C MARYSIA
c endif !sscale
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
& +a33*muij(4)
c write (iout,*) 'i',i,' j',j,itype(i),itype(j),
c & ' eel_loc_ij',eel_loc_ij
-c write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
+C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
C Calculate patrial derivative for theta angle
#ifdef NEWCORR
geel_loc_ij=a22*gmuij1(1)
etheta=0.0D0
c write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
+c write (iout,*) "ebend: i=",i
+c call flush_(iout)
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
+c write (iout,*) "Exit loop"
+c call flush_(iout)
ethetacnstr=0.0d0
-C print *,ithetaconstr_start,ithetaconstr_end,"TU"
- do i=ithetaconstr_start,ithetaconstr_end
+c write (iout,*) ithetaconstr_start,ithetaconstr_end,"TU"
+c call flush_(iout)
+ do i=max0(ithetaconstr_start,1),ithetaconstr_end
itheta=itheta_constr(i)
thetiii=theta(itheta)
difi=pinorm(thetiii-theta_constr0(i))
& gloc(itheta+nphi-2,icg)
endif
enddo
+c write (iout,*) "Exit ebend"
+c call flush_(iout)
C Ufff.... We've done all this!!!
return
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+ if (i.le.2) cycle
c print *,i,itype(i-1),itype(i),itype(i-2)
if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
& .or.itype(i).eq.ntyp1) cycle
sinkt(k)=dsin(k*theti2)
enddo
C print *,ethetai
+ if (i.eq.3) then
+ phii=0.0d0
+ ityp1=nthetyp+1
+ do k=1,nsingle
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ else
+
if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
sinph1(k)=0.0d0
enddo
endif
+ endif
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
C now constrains
ethetacnstr=0.0d0
C print *,ithetaconstr_start,ithetaconstr_end,"TU"
- do i=ithetaconstr_start,ithetaconstr_end
+ do i=max0(ithetaconstr_start,1),ithetaconstr_end
itheta=itheta_constr(i)
thetiii=theta(itheta)
difi=pinorm(thetiii-theta_constr0(i))
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
-c lprn=.true.
+c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
C ANY TWO ARE DUMMY ATOMS in row CYCLE
if (lprn) then
write (iout,'(a)') 'Contact function values before RECEIVE:'
do i=nnt,nct-2
- write (iout,'(2i3,50(1x,i2,f5.2))')
+ write (iout,'(2i3,50(1x,i3,f5.2))')
& i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
& j=1,num_cont_hb(i))
enddo
jp1=iabs(j1)
c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
c & ' jj=',jj,' kk=',kk
+c call flush(iout)
if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0
& .or. j.lt.0 .and. j1.gt.0) .and.
& (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then
enddo ! kk
do kk=1,num_conti
j1=jcont_hb(kk,i)
-c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
-c & ' jj=',jj,' kk=',kk
+c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+c & ' jj=',jj,' kk=',kk
+c call flush(iout)
if (j1.eq.j+1) then
C Contacts I-J and (I+1)-J occur simultaneously.
C The system loses extra energy.
include 'COMMON.CONTACTS'
include 'COMMON.CHAIN'
include 'COMMON.CONTROL'
+ include 'COMMON.TORSION'
double precision gx(3),gx1(3)
integer num_cont_hb_old(maxres)
logical lprn,ldone
C Set lprn=.true. for debugging
lprn=.false.
eturn6=0.0d0
+c write (iout,*) "MULTIBODY_EELLO"
+c call flush(iout)
#ifdef MPI
do i=1,nres
num_cont_hb_old(i)=num_cont_hb(i)
if (lprn) then
write (iout,'(a)') 'Contact function values before RECEIVE:'
do i=nnt,nct-2
- write (iout,'(2i3,50(1x,i2,f5.2))')
+ write (iout,'(2i3,50(1x,i3,f5.2))')
& i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
& j=1,num_cont_hb(i))
enddo
+ call flush(iout)
endif
- call flush(iout)
do i=1,ntask_cont_from
ncont_recv(i)=0
enddo
if (lprn) then
write (iout,'(a)') 'Contact function values:'
do i=nnt,nct-2
- write (iout,'(2i3,50(1x,i2,5f6.3))')
+ write (iout,'(2i3,50(1x,i3,5f6.3))')
& i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i),
& ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i))
enddo
+ write (iout,*) "itortyp"
+ do i=1,nres
+ write (iout,*) i,itype(i),itortyp(itype(i))
+ enddo
+ call flush(iout)
endif
ecorr=0.0D0
ecorr5=0.0d0
do kk=1,num_conti1
j1=jcont_hb(kk,i1)
jp1=iabs(j1)
-c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
-c & ' jj=',jj,' kk=',kk
+ if (lprn) then
+ write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+ & ' jj=',jj,' kk=',kk
+ call flush(iout)
+ endif
c if (j1.eq.j+1 .or. j1.eq.j-1) then
if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0
& .or. j.lt.0 .and. j1.gt.0) .and.
c write (iout,'(i5,3f10.5)')
c & iii,(gradcorr5(jjj,iii),jjj=1,3)
c enddo
+c write (iout,*) "ecorr4"
+c call flush(iout)
+c write (iout,*) "eello5:",i,jp,i+1,jp1,jj,kk,
+c & itype(jp),itype(i+1),itype(jp1),
+c & itortyp(itype(jp)),itortyp(itype(i+1)),itortyp(itype(jp1))
+c call flush(iout)
if (wcorr5.gt.0.0d0)
& ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk)
c write (iout,*) "gradcorr5 after eello5"
2 'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk)
cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6
cd write(2,*)'ijkl',i,jp,i+1,jp1
+c write (iout,*) "ecorr5"
+c call flush(iout)
if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3
& .or. wturn6.eq.0.0d0))then
cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk)
if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
1 'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk)
+c write (iout,*) "ecorr6"
+c call flush(iout)
cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
cd & 'ecorr6=',ecorr6
cd write (iout,'(4e15.5)') sred_geom,
if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
1 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk)
cd write (2,*) 'multibody_eello:eturn6',eturn6
+c write (iout,*) "ecorr4"
+c call flush(iout)
endif
ENDIF
1111 continue
endif
+ if (energy_dec) call flush(iout)
enddo ! kk
enddo ! jj
enddo ! 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))
+c itk=itortyp(itype(k))
+c itl=itortyp(itype(l))
+c itj=itortyp(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))
+c itk=itortyp(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
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
C energy moment and not to the cluster cumulant.
- iti=itortyp(itype(i))
- if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
- else
- itj1=ntortyp
- endif
- itk=itortyp(itype(k))
- itk1=itortyp(itype(k+1))
- if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
- else
- itl1=ntortyp
- endif
+c iti=itortyp(itype(i))
+c if (j.lt.nres-1) then
+c itj1=itortyp(itype(j+1))
+c else
+c itj1=ntortyp
+c endif
+c itk=itortyp(itype(k))
+c itk1=itortyp(itype(k+1))
+c if (l.lt.nres-1) then
+c itl1=itortyp(itype(l+1))
+c else
+c itl1=ntortyp
+c endif
#ifdef MOMENT
s1=dip(4,jj,i)*dip(4,kk,k)
#endif
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))
- if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
- else
- itj1=ntortyp
- endif
- itk=itortyp(itype(k))
- if (k.lt.nres-1) then
- itk1=itortyp(itype(k+1))
- else
- itk1=ntortyp
- endif
- itl=itortyp(itype(l))
- if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
- else
- itl1=ntortyp
- endif
+c iti=itortyp(itype(i))
+c itj=itortyp(itype(j))
+c if (j.lt.nres-1) then
+c itj1=itortyp(itype(j+1))
+c else
+c itj1=ntortyp
+c endif
+c itk=itortyp(itype(k))
+c if (k.lt.nres-1) then
+c itk1=itortyp(itype(k+1))
+c else
+c itk1=ntortyp
+c endif
+c itl=itortyp(itype(l))
+c if (l.lt.nres-1) then
+c itl1=itortyp(itype(l+1))
+c else
+c itl1=ntortyp
+c 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,
cd & ' itl',itl,' itl1',itl1
include 'COMMON.INTERACT'
C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
double precision div77_81/0.974996043d0/,
- &div4_81/0.2222222222d0/
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
C the vector between center of side_chain and peptide group
double precision pep_side(3),long,side_calf(3),
- &pept_group(3)
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
C the line belowe needs to be changed for FGPROC>1
do i=1,nres-1
if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
C the line below has to be changed for FGPROC>1
VolumeTotal=0.0
do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
dist_pep_side=0.0
dist_side_calf=0.0
do j=1,3
C first lets set vector conecting the ithe side-chain with kth side-chain
- pep_side(j)=c(k+nres,j)-(c(i,j)+c(i+1,j))/2.0d0
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
C and vector conecting the side-chain with its proper calfa
- side_calf(j)=c(k+nres,j)-c(k,j)
- pept_group(j)=c(i,j)-c(i+1,j)
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
C lets have their lenght
dist_pep_side=pep_side(j)**2+dist_pep_side
dist_side_calf=dist_side_calf+side_calf(j)**2
enddo
dist_pep_side=dsqrt(dist_pep_side)
dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
C now sscale fraction
sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
C now sscale
if (sh_frac_dist.le.0.0) cycle
C If we reach here it means that this side chain reaches the shielding sphere
ishield_list(i)=ishield_list(i)+1
C ishield_list is a list of non 0 side-chain that contribute to factor gradient
C this list is essential otherwise problem would be O3
- shield_list(ishield_list)=k
+ shield_list(ishield_list(i),i)=k
C Lets have the sscale value
if (sh_frac_dist.gt.1.0) then
scale_fac_dist=1.0d0
else
scale_fac_dist=-sh_frac_dist*sh_frac_dist
& *(2.0*sh_frac_dist-3.0d0)
- fac_help_scale=6.0*(scale_fac_dist-scale_fac_dist**2)
+ fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
& /dist_pep_side/buff_shield*0.5
C remember for the final gradient multiply sh_frac_dist_grad(j)
C for side_chain by factor -2 !
do j=1,3
sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
enddo
endif
+C if ((i.eq.3).and.(k.eq.2)) then
+C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C & ,"TU"
+C endif
+
C this is what is now we have the distance scaling now volume...
short=short_r_sidechain(itype(k))
long=long_r_sidechain(itype(k))
- costhet=1.0d0/dsqrt(1+short**2/dist_pep_side**2)
+ costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
C now costhet_grad
- costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side
+C costhet=0.0d0
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+C costhet_fac=0.0d0
do j=1,3
costhet_grad(j)=costhet_fac*pep_side(j)
enddo
do j=1,3
pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
enddo
- fac_alfa_sin=1.0-(pep_side0pept_group/
- & (dist_pep_side*dist_side_calf))**2
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0-cosalfa**2
fac_alfa_sin=dsqrt(fac_alfa_sin)
rkprim=fac_alfa_sin*(long-short)+short
- cosphi=1.0d0/dsqrt(1+rkprim**2/dist_pep_side**2)
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+
+ cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+ enddo
+
VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
& /VSolvSphere_div
+C now the gradient...
+C grad_shield is gradient of Calfa for peptide groups
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)
+C gradient po costhet
+ &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+ &-scale_fac_dist*(cosphi_grad_long(j))
+ &/(1.0-cosphi) )*div77_81
+ &*VofOverlap
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*(-2.0d0)
+ & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+ & +scale_fac_dist*(cosphi_grad_long(j))
+ & *2.0d0/(1.0-cosphi))
+ & *div77_81*VofOverlap
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*cosphi_grad_loc(j)
+ & *2.0d0/(1.0-cosphi)
+ & *div77_81*VofOverlap
+ enddo
VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
-C if ((cosphi.le.0.0).or.(costhet.le.0.0)) write(iout,*) "ERROR",
-C & cosphi,costhet
-C now should be fac_side_grad(k) which will be gradient of factor k which also
-C affect the gradient of peptide group i fac_pept_grad(i) and i+1
- write(2,*) "myvolume",VofOverlap,VSolvSphere_div,VolumeTotal
- enddo
-C write(2,*) "TOTAL VOLUME",i,VolumeTotal
-C the scaling factor of the shielding effect
+ enddo
fac_shield(i)=VolumeTotal*div77_81+div4_81
- write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
enddo
return
end
projects / unres.git / blobdiff