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 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
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
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
&*VofOverlap
C grad_shield_side is Cbeta sidechain gradient
grad_shield_side(j,ishield_list(i),i)=
- & (sh_frac_dist_grad(j)*-2.0d0
+ & (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))
projects / unres.git / blobdiff