projects
/
unres.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
update
[unres.git]
/
source
/
unres
/
src_MD
/
energy_p_new_barrier.F
diff --git
a/source/unres/src_MD/energy_p_new_barrier.F
b/source/unres/src_MD/energy_p_new_barrier.F
index
07445bc
..
4bccdf3
100644
(file)
--- a/
source/unres/src_MD/energy_p_new_barrier.F
+++ b/
source/unres/src_MD/energy_p_new_barrier.F
@@
-131,6
+131,11
@@
C
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+cmc
+cmc Sep-06: egb takes care of dynamic ss bonds too
+cmc
+c if (dyn_ss) call dyn_set_nss
+
c print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
#ifdef MPI
c print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
#ifdef MPI
@@
-326,6
+331,7
@@
C
energia(23)=evdw_m
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
energia(23)=evdw_m
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
+ if (dyn_ss) call dyn_set_nss
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
#ifdef MPI
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
#ifdef MPI
@@
-1552,6
+1558,7
@@
C
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
+ include 'COMMON.SBRIDGE'
logical lprn
evdw=0.0D0
ccccc energy_dec=.false.
logical lprn
evdw=0.0D0
ccccc energy_dec=.false.
@@
-1580,6
+1587,12
@@
C Calculate SC interaction energy.
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,' ss'
+ ELSE
ind=ind+1
itypj=itype(j)
c dscj_inv=dsc_inv(itypj)
ind=ind+1
itypj=itype(j)
c dscj_inv=dsc_inv(itypj)
@@
-1666,9
+1679,10
@@
c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
& evdwij
endif
& evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw',i,j,evdwij
-
+ if (energy_dec) then
+ write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+ call flush(iout)
+ endif
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
@@
-1689,6
+1703,7
@@
C Calculate angular part of the gradient.
#else
call sc_grad
#endif
#else
call sc_grad
#endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
enddo ! j
enddo ! iint
enddo ! i
@@
-3009,6
+3024,9
@@
C
C Loop over i,i+2 and i,i+3 pairs of the peptide groups
C
do i=iturn3_start,iturn3_end
C Loop over i,i+2 and i,i+3 pairs of the peptide groups
C
do i=iturn3_start,iturn3_end
+C if (itype(i).eq.21 .or. itype(i+1).eq.21
+C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21)
+C & cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
@@
-3024,6
+3042,10
@@
C
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
+C if (itype(i).eq.21 .or. itype(i+1).eq.21
+C & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21
+C & .or. itype(i+5).eq.21)
+C & cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
@@
-3042,6
+3064,8
@@
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
+C if (itype(i).eq.21 .or. itype(i+1).eq.21
+C &.or.itype(i+2)) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
@@
-3054,6
+3078,8
@@
c
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
+C if (itype(j).eq.21 .or. itype(j+1).eq.21
+C &.or.itype(j+2)) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
@@
-4266,9
+4292,15
@@
c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
+ if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
+ if (ii.gt.nres
+ & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+ endif
cd write (iout,*) "eij",eij
else if (ii.gt.nres .and. jj.gt.nres) then
c Restraints from contact prediction
cd write (iout,*) "eij",eij
else if (ii.gt.nres .and. jj.gt.nres) then
c Restraints from contact prediction
@@
-4462,6
+4494,8
@@
c
do i=ibondp_start,ibondp_end
diff = vbld(i)-vbldp0
c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
do i=ibondp_start,ibondp_end
diff = vbld(i)-vbldp0
c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+ if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+ & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
@@
-4480,6
+4514,12
@@
c
diff=vbld(i+nres)-vbldsc0(1,iti)
c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
c & AKSC(1,iti),AKSC(1,iti)*diff*diff
diff=vbld(i+nres)-vbldsc0(1,iti)
c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
c & AKSC(1,iti),AKSC(1,iti)*diff*diff
+ if (energy_dec) then
+ write (iout,*)
+ & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+ & AKSC(1,iti),AKSC(1,iti)*diff*diff
+ call flush(iout)
+ endif
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@
-4758,6
+4798,8
@@
C
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
logical lprn /.false./, lprn1 /.false./
etheta=0.0D0
do i=ithet_start,ithet_end
+ if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+ &(itype(i).eq.ntyp1)) cycle
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
@@
-4767,7
+4809,8
@@
C
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3) then
+C if (i.gt.3) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
@@
-4781,13
+4824,13
@@
C
enddo
else
phii=0.0d0
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres) then
+ if ((i.lt.nres).and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
@@
-4802,7
+4845,7
@@
C
enddo
else
phii1=0.0d0
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp(itype(i))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
@@
-4908,10
+4951,12
@@
C
enddo
enddo
10 continue
enddo
enddo
10 continue
- if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
- & i,theta(i)*rad2deg,phii*rad2deg,
+ if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)')
+ & 'ebe', i,theta(i)*rad2deg,phii*rad2deg,
& phii1*rad2deg,ethetai
etheta=etheta+ethetai
& phii1*rad2deg,ethetai
etheta=etheta+ethetai
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ & 'ebend',i,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=wang*dethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=wang*dethetai
@@
-5853,12
+5898,14
@@
C 6/23/01 Compute double torsional energy
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.TORCNSTR'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.TORCNSTR'
+ include 'COMMON.CONTROL'
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
etors_d=0.0D0
do i=iphid_start,iphid_end
logical lprn
C Set lprn=.true. for debugging
lprn=.false.
c lprn=.true.
etors_d=0.0D0
do i=iphid_start,iphid_end
+ etors_d_ii=0.0D0
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
itori2=itortyp(itype(i))
@@
-5877,6
+5924,8
@@
c lprn=.true.
sinphi2=dsin(j*phii1)
etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
& v2cij*cosphi2+v2sij*sinphi2
sinphi2=dsin(j*phii1)
etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
& v2cij*cosphi2+v2sij*sinphi2
+ if (energy_dec) etors_d_ii=etors_d_ii+
+ & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
@@
-5892,12
+5941,17
@@
c lprn=.true.
sinphi1m2=dsin(l*phii-(k-l)*phii1)
etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
& v1sdij*sinphi1p2+v2sdij*sinphi1m2
sinphi1m2=dsin(l*phii-(k-l)*phii1)
etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
& v1sdij*sinphi1p2+v2sdij*sinphi1m2
+ if (energy_dec) etors_d_ii=etors_d_ii+
+ & v1cdij*cosphi1p2+v2cdij*cosphi1m2+
+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2
gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
& -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
& -v1cdij*sinphi1p2+v2cdij*sinphi1m2)
enddo
enddo
gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
& -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
& -v1cdij*sinphi1p2+v2cdij*sinphi1m2)
enddo
enddo
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ & 'etor_d',i,etors_d_ii
gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
c write (iout,*) "gloci", gloc(i-3,icg)
gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
c write (iout,*) "gloci", gloc(i-3,icg)
@@
-5934,18
+5988,22
@@
c lprn=.true.
c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
esccor=0.0D0
do i=itau_start,itau_end
c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
esccor=0.0D0
do i=itau_start,itau_end
- esccor_ii=0.0D0
+C do i=42,42
+
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
phii=phi(i)
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
phii=phi(i)
+
cccc Added 9 May 2012
cc Tauangle is torsional engle depending on the value of first digit
c(see comment below)
cc Omicron is flat angle depending on the value of first digit
c(see comment below)
cccc Added 9 May 2012
cc Tauangle is torsional engle depending on the value of first digit
c(see comment below)
cc Omicron is flat angle depending on the value of first digit
c(see comment below)
-
+C print *,i,tauangle(1,i)
do intertyp=1,3 !intertyp
do intertyp=1,3 !intertyp
+ esccor_ii=0.0D0
cc Added 09 May 2012 (Adasko)
cc Intertyp means interaction type of backbone mainchain correlation:
c 1 = SC...Ca...Ca...Ca
cc Added 09 May 2012 (Adasko)
cc Intertyp means interaction type of backbone mainchain correlation:
c 1 = SC...Ca...Ca...Ca
@@
-5967,9
+6025,13
@@
c 3 = SC...Ca...Ca...SCi
v2ij=v2sccor(j,intertyp,isccori,isccori1)
cosphi=dcos(j*tauangle(intertyp,i))
sinphi=dsin(j*tauangle(intertyp,i))
v2ij=v2sccor(j,intertyp,isccori,isccori1)
cosphi=dcos(j*tauangle(intertyp,i))
sinphi=dsin(j*tauangle(intertyp,i))
+ if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+ if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)')
+ & 'esccor',i,intertyp,esccor_ii
+C print *,i,tauangle(1,i),gloci
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
c &gloc_sc(intertyp,i-3,icg)
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
c &gloc_sc(intertyp,i-3,icg)
@@
-5982,7
+6044,8
@@
c &gloc_sc(intertyp,i-3,icg)
enddo !intertyp
enddo
c do i=1,nres
enddo !intertyp
enddo
c do i=1,nres
-c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg)
+c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc_sc(2,i,icg),
+c & gloc_sc(3,i,icg)
c enddo
return
end
c enddo
return
end