integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc,
& nres0,nstart_seq,nchain,chain_length,chain_border,iprzes,
- & chain_border1,ireschain,tabpermchain,npermchain,afmend,afmbeg,
- & nres_chomo,nmodel_start
+ & chain_border1,ireschain,tabpermchain,npermchain,nequiv,
+ & nchain_group,iequiv,mapchain,afmend,afmbeg,
+ & nres_chomo,nmodel_start,nran_start
double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r,
& prod,rt,dc_work,cref,crefjlee,dc_norm2,velAFMconst,
& totTafm,chomo
common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
& xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
& dc_norm2(3,0:maxres2),
- & dc_work(MAXRES6),nres,nres0
+ & dc_work(MAXRES6),nres,nres0,nran_start
common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres),
& rt(3,3,maxres)
common /refstruct/ cref(3,maxres2+2),
& crefjlee(3,maxres2+2),
& nsup,nstart_sup,nstart_seq,iprzes,
& chain_length(maxchain),npermchain,ireschain(maxres),
- & tabpermchain(maxchain,maxperm),
+ & tabpermchain(maxchain,maxperm),nchain_group,
+ & iequiv(maxchain,maxres),nequiv(maxchain),mapchain(maxchain),
& chain_border(2,maxchain),chain_border1(2,maxchain),nchain
common /from_zscore/ nz_start,nz_end,iz_sc
double precision boxxsize,boxysize,boxzsize,enecut,sscut,
! and output level.
!... energy_dec = .true. means print energy decomposition matrix
integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
- & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
+ & inprint,i2ndstr,mucadyn,constr_dist,symetr,npermut,
+ & AFMlog,selfguide,
& shield_mode,tor_mode,tubelog,constr_homology,homol_nset,
& iprint
!... minim = .true. means DO minimization.
& unres_pdb,out_cart,out_int,vdisulf,searchsc,lmuca,dccart,mucadyn,
& extconf,out1file,gmatout,selfguide,AFMlog,shield_mode,tor_mode,
& tubelog,constr_dist,gnorm_check,gradout,split_ene,
- & with_theta_constr,with_dihed_constr,symetr,usampl,loc_qlike,
+ & with_theta_constr,with_dihed_constr,symetr,npermut,usampl,loc_
+ & qlike,
& adaptive,constr_homology,homol_nset,read2sigma,start_from_model,
& read_homol_frag,out_template_coord,out_template_restr
& iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,iatel_e_vdw,
& iatscp_s,iatscp_e,ispp,iscp
C 3/26/20 Interaction lists
+C 10/26/20 Upgraded to interaction lists for RESPA
integer newcontlisti(maxint_res*maxres),
+ & newcontlisti_long(maxint_res*maxres),
+ & newcontlisti_short(maxint_res*maxres),
& newcontlistj(maxint_res*maxres),
+ & newcontlistj_long(maxint_res*maxres),
+ & newcontlistj_short(maxint_res*maxres),
& newcontlistppi(maxint_res*maxres),
& newcontlistppj(maxint_res*maxres),
- & newcontlistpp_vdwi(maxint_res*maxres),
- & newcontlistpp_vdwj(maxint_res*maxres),
+ & newcontlistpp_vdwi_short(maxint_res*maxres),
+ & newcontlistpp_vdwj_short(maxint_res*maxres),
& newcontlistscpi(2*maxint_res*maxres),
+ & newcontlistscpi_long(2*maxint_res*maxres),
+ & newcontlistscpi_short(2*maxint_res*maxres),
& newcontlistscpj(2*maxint_res*maxres),
+ & newcontlistscpj_long(2*maxint_res*maxres),
+ & newcontlistscpj_short(2*maxint_res*maxres),
& g_listscsc_start,g_listscsc_end,g_listpp_start,g_listpp_end,
- & g_listpp_vdw_start,g_listpp_vdw_end,g_listscp_start,g_listscp_end
- common /interact_list/newcontlisti,newcontlistj,g_listscsc_start,
- & g_listscsc_end,newcontlistppi,newcontlistppj,g_listpp_start,
- & g_listpp_end,newcontlistpp_vdwi,newcontlistpp_vdwj,
- & g_listpp_vdw_start,g_listpp_vdw_end,
- & newcontlistscpi,newcontlistscpj,g_listscp_start,
- & g_listscp_end
+ & g_listscp_start,g_listscp_end,
+ & g_listscsc_start_long,g_listscsc_end_long,
+ & g_listscp_start_long,g_listscp_end_long,
+ & g_listscsc_start_short,g_listscsc_end_short,
+ & g_listpp_vdw_start_short,g_listpp_vdw_end_short,
+ & g_listscp_start_short,g_listscp_end_short
+ common /interact_list/
+ & newcontlisti,newcontlisti_long,newcontlisti_short,
+ & newcontlistj,newcontlistj_long,newcontlistj_short,
+ & g_listscsc_start,g_listscsc_start_long,g_listscsc_start_short,
+ & g_listscsc_end,g_listscsc_end_long,g_listscsc_end_short,
+ & newcontlistppi,newcontlistppj,g_listpp_start,
+ & g_listpp_end,
+ & newcontlistpp_vdwi_short,
+ & newcontlistpp_vdwj_short,
+ & g_listpp_vdw_start_short,
+ & g_listpp_vdw_end_short,
+ & newcontlistscpi,newcontlistscpi_long,newcontlistscpi_short,
+ & newcontlistscpj,newcontlistscpj_long,newcontlistscpj_short,
+ & g_listscp_start,g_listscp_start_long,g_listscp_start_short,
+ & g_listscp_end,g_listscp_end_long,g_listscp_end_short
C 12/1/95 Array EPS included in the COMMON block.
double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,alp,
& sigma0,sigii,
& dvmax,damax,edriftmax
integer n_timestep,ntwx,ntwe,lang,count_reset_moment,
& count_reset_vel,ntime_split,ntime_split0,
- & maxtime_split,itime_mat,imatupdate
+ & maxtime_split,itime_mat,imatupdate,irest_freq
logical large,print_compon,tbf,rest,reset_moment,reset_vel,
& rattle,mdpdb,RESPA,preminim
common /mdpar/ v_ini,d_time,d_time0,t_bath,
& tau_bath,dvmax,damax,n_timestep,mdpdb,
& ntime_split,ntime_split0,maxtime_split,itime_mat,imatupdate,
- & ntwx,ntwe,lang,large,print_compon,tbf,rest,preminim,
+ & ntwx,ntwe,irest_freq,lang,large,print_compon,tbf,rest,preminim,
& reset_moment,reset_vel,count_reset_moment,count_reset_vel,
& rattle,RESPA
! Basic quantities
& time_vec,time_mat,time_ginvmult,time_fricmatmult,time_fric,
& time_scatter_fmat,time_scatter_ginv,
& time_fsample,time_scatter_fmatmult,time_scatter_ginvmult,
- & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS
+ & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS,time_list,
+ & time_evdw,time_eelec,time_escp,time_evdw_short,time_eelec_short,
+ & time_escp_short,time_evdw_long,time_eelec_long,time_escp_long,
+ & time_escpcalc,time_escpsetup
common /timing/ t_init,t_MDsetup,t_langsetup,
& t_MD,t_enegrad,t_sdsetup,time_bcast,time_reduce,time_gather,
& time_sendrecv,time_scatter,time_barrier_e,time_barrier_g,
& time_vec,time_mat,time_ginvmult,time_fricmatmult,time_fric,
& time_fsample,time_scatter_fmatmult,time_scatter_ginvmult,
& time_scatter_fmat,time_scatter_ginv,
- & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS
+ & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS,time_list,
+ & time_evdw,time_eelec,time_escp,time_evdw_short,time_eelec_short,
+ & time_escp_short,time_evdw_long,time_eelec_long,time_escp_long,
+ & time_escpcalc,time_escpsetup
- integer ndih_constr,idih_constr(maxdih_constr),ntheta_constr,
+ integer ndih_constr,idih_constr(maxdih_constr),
+ & iconstr_dih(maxres),ntheta_constr,
& itheta_constr(maxdih_constr)
integer ndih_nconstr,idih_nconstr(maxdih_constr)
integer idihconstr_start,idihconstr_end,ithetaconstr_start,
& vpsipred(3,maxdih_constr),sdihed(2,maxdih_constr),wdihc
common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange,
& for_thet_constr,vpsipred,sdihed,wdihc,
- & ndih_constr,idih_constr,
+ & ndih_constr,idih_constr,iconstr_dih,
& ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end,
& ntheta_constr,itheta_constr,ithetaconstr_start,
& ithetaconstr_end,raw_psipred
parameter (max_cg_procs=maxprocs)
C Max. number of AA residues
integer maxres
- parameter (maxres=20000)
+c parameter (maxres=70000)
+ parameter (maxres=1000)
C Max. number of AA residues per chain
integer maxres_chain
parameter (maxres_chain=1200)
& mmaxres2_chain=maxres2_chain*(maxres2_chain+1)/2)
C Max number of symetric chains
integer maxchain
- parameter (maxchain=50)
+ parameter (maxchain=250)
integer maxperm
- parameter (maxperm=5040)
+c parameter (maxperm=5040)
+ parameter (maxperm=3628800)
C Max. number of variables
integer maxvar
parameter (maxvar=6*maxres)
C Max. number of derivatives of virtual-bond and side-chain vectors in theta
C or phi.
integer maxdim
-c parameter (maxdim=(maxres_chain-1)*(maxres_chain-2)/2)
- parameter (maxdim=(maxres-1)*(maxres-2)/2)
+ parameter (maxdim=(maxres_chain-1)*(maxres_chain-2)/2)
+c parameter (maxdim=(maxres-1)*(maxres-2)/2)
C Max. number of SC contacts
integer maxcont
parameter (maxcont=12*maxres)
parameter (maxcont_res=200)
C Max. number of distance/contact-distance restraints
integer maxdim_cont
-c parameter (maxdim_cont=maxres*maxcont_res)
- parameter (maxdim_cont=maxres*1000)
+ parameter (maxdim_cont=maxres*maxcont_res)
C Number of AA types (at present only natural AA's will be handled
integer ntyp,ntyp1
parameter (ntyp=24,ntyp1=ntyp+1)
t_enegrad=0.0d0
t_sdsetup=0.0d0
write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
+ write (iout,*) "Restart frequency:",irest_freq
#ifdef MPI
tt0=MPI_Wtime()
#else
call cartout(totT)
endif
endif
- if (rstcount.eq.1000.or.itime.eq.n_timestep) then
+ if (rstcount.eq.irest_freq.or.itime.eq.n_timestep) then
open(irest2,file=rest2name,status='unknown')
write(irest2,*) totT,EK,potE,totE,t_bath
do i=0,2*nres-1
& ' End of MD calculation '
#ifdef TIMING_ENE
write (iout,*) "time for etotal",t_etotal," elong",t_elong,
- & " eshort",t_eshort
+ & " eshort",t_eshort," list",time_list
write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
& " time_fricmatmult",time_fricmatmult," time_fsample ",
& time_fsample
call zerograd
call etotal(potEcomp)
! AL 4/17/17: Reduce the steps if NaNs occurred.
- if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
+ if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0))) then
d_time=d_time/2
cycle
endif
double precision etot
logical fail
integer i_start_models(0:nodes-1)
+ double precision potEcomp_long(0:Max_Ene)
write (iout,*) "init_MD INDPDB",indpdb
d_time0=d_time
c write(iout,*) "d_time", d_time
if(me.eq.king.or..not.out1file)then
write (iout,*) "Initial velocities"
do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+ write (iout,'(i7,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
& (d_t(j,i+nres),j=1,3)
enddo
endif
call zerograd
call etotal(potEcomp)
if (large) call enerprint(potEcomp)
+ if (RESPA) call etotal_long(potEcomp_long)
#ifdef TIMING_ENE
#ifdef MPI
t_etotal=t_etotal+MPI_Wtime()-tt0
call cartgrad
call lagrangian
call max_accel
- if (amax*d_time .gt. dvmax) then
+ if (amax*d_time .gt. dvmax .and. .not. respa) then
d_time=d_time*dvmax/amax
if(me.eq.king.or..not.out1file) write (iout,*)
& "Time step reduced to",d_time,
& ' centroid coordinates'/
1 ' ', 7X,'X',11X,'Y',11X,'Z',
& 10X,'X',11X,'Y',11X,'Z')
- 110 format (a,'(',i4,')',6f12.5)
+ 110 format (a,'(',i5,')',6f12.5)
return
end
subroutine chain_symmetry(nchain,nres,itype,chain_border,
- & chain_length,npermchain,tabpermchain)
+ & chain_length,npermchain,tabpermchain,nchain_group,nequiv,
+ & iequiv,chaingroup)
c
c Determine chain symmetry. nperm is the number of permutations and
c tabperchain contains the allowed permutations of the chains.
& chain_length(nchain),itemp(maxchain),
& npermchain,tabpermchain(maxchain,maxperm),
& tabperm(maxchain,maxperm),mapchain(maxchain),
- & iequiv(maxchain,maxres),iflag(maxres)
+ & chaingroup(maxchain),iequiv(maxchain,maxres),iflag(maxres)
integer i,j,k,l,ii,nchain_group,nequiv(maxchain),iieq,
& nperm,npermc,ind
+ logical lprn /.false./
if (nchain.eq.1) then
npermchain=1
tabpermchain(1,1)=1
c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
+ nchain_group=1
+ iequiv(1,1)=1
+ chaingroup(1)=1
return
endif
c
do j=1,nequiv(i)
ind=ind+1
mapchain(ind)=iequiv(j,i)
+ chaingroup(ind)=i
enddo
enddo
write (iout,*) "mapchain"
do i=1,nchain
write (iout,*) i,mapchain(i)
enddo
+ write (iout,*) "chaingroup"
+ do i=1,nchain
+ write (iout,*) i,chaingroup(i)
+ enddo
+ if (npermut.eq.0) then
+ npermchain=1
+ do i=1,nchain
+ tabpermchain(i,1)=i
+ enddo
+ return
+ endif
ii=0
do i=1,nchain_group
call permut(nequiv(i),nperm,tabperm)
enddo
enddo
write(iout,*) "Number of chain permutations",npermchain
+ if (lprn) then
write(iout,*) "Permutations"
do i=1,npermchain
write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
enddo
+ endif
return
end
c---------------------------------------------------------------------
include 'COMMON.VAR'
cost=dcos(theta(3))
sint=dsin(theta(3))
- t(1,1,1)=-cost
+ t(1,1,1)=cost
t(1,2,1)=-sint
t(1,3,1)= 0.0D0
- t(2,1,1)=-sint
- t(2,2,1)= cost
+ t(2,1,1)=sint
+ t(2,2,1)=cost
t(2,3,1)= 0.0D0
t(3,1,1)= 0.0D0
t(3,2,1)= 0.0D0
C Locate CA(i) and SC(i-1)
C
implicit none
- integer i,j
+ integer i,j,k
double precision theti,phii,cost,sint,cosphi,sinphi
include 'DIMENSIONS'
include 'COMMON.CHAIN'
sint=dsin(theti)
cosphi=dcos(phii)
sinphi=dsin(phii)
+c write (iout,*) "locate_next_res i",i
* Define the matrices of the rotation about the virtual-bond valence angles
* theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
* program), R(i,j,k), and, the cumulative matrices of rotation RT
t(1,1,i-2)=-cost
- t(1,2,i-2)=-sint
+ t(1,2,i-2)= sint
t(1,3,i-2)= 0.0D0
t(2,1,i-2)=-sint
- t(2,2,i-2)= cost
+ t(2,2,i-2)=-cost
t(2,3,i-2)= 0.0D0
t(3,1,i-2)= 0.0D0
t(3,2,i-2)= 0.0D0
r(1,2,i-2)= 0.0D0
r(1,3,i-2)= 0.0D0
r(2,1,i-2)= 0.0D0
- r(2,2,i-2)=-cosphi
- r(2,3,i-2)= sinphi
+ r(2,2,i-2)= cosphi
+ r(2,3,i-2)=-sinphi
r(3,1,i-2)= 0.0D0
r(3,2,i-2)= sinphi
r(3,3,i-2)= cosphi
rt(1,1,i-2)=-cost
- rt(1,2,i-2)=-sint
+ rt(1,2,i-2)= sint
rt(1,3,i-2)=0.0D0
- rt(2,1,i-2)=sint*cosphi
+ rt(2,1,i-2)=-sint*cosphi
rt(2,2,i-2)=-cost*cosphi
- rt(2,3,i-2)=sinphi
+ rt(2,3,i-2)=-sinphi
rt(3,1,i-2)=-sint*sinphi
- rt(3,2,i-2)=cost*sinphi
- rt(3,3,i-2)=cosphi
+ rt(3,2,i-2)=-cost*sinphi
+ rt(3,3,i-2)= cosphi
call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
+c write (iout,*) "prod",i-2
+c do j=1,3
+c write (iout,*) (prod(j,k,i-2),k=1,3)
+c enddo
+c write (iout,*) "prod",i-1
+c do j=1,3
+c write (iout,*) (prod(j,k,i-1),k=1,3)
+c enddo
do j=1,3
dc_norm(j,i-1)=prod(j,1,i-1)
dc(j,i-1)=vbld(i)*prod(j,1,i-1)
c(j,i)=c(j,i-1)+dc(j,i-1)
enddo
+c write (iout,*) "dc",i-1,(dc(j,i-1),j=1,3)
+c write (iout,*) "c",i,(dc(j,i),j=1,3)
cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
C
C Now calculate the coordinates of SC(i-1)
double precision dnorm1,dnorm2,be
double precision time00
double precision dist,alpha,beta
+ double precision time01
if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
#ifdef TIMING
time01=MPI_Wtime()
i2=icont(2,i)
it1=itype(i1)
it2=itype(i2)
- write (iout,'(i5,2x,a,i5,2x,a,i4)')
+ write (iout,'(i10,2x,a,i7,2x,a,i7)')
& i,restyp(it1),i1,restyp(it2),i2
enddo
endif
ii1=iharp(3,i)
jj1=iharp(4,i)
write (iout,*)
- write (iout,'(20(a,i5,1x))') (restyp(itype(k)),k,k=i1,ii1)
- write (iout,'(20(a,i5,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
+ write (iout,'(20(a,i7,1x))') (restyp(itype(k)),k,k=i1,ii1)
+ write (iout,'(20(a,i7,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
c do k=jj1,j1,-1
c write (iout,'(a,i3,$)') restyp(itype(k)),k
c enddo
edfadis=0
gdfad=0.0d0
-c write (2,*) "edfad",idfadis_start,idfadis_end
do i=idfadis_start,idfadis_end
iatm1=idislis(1,i)+ishiftca
include 'COMMON.NAMES'
logical lprint
double precision elpp_6(2,2),elpp_3(2,2),ael6_(2,2),ael3_(2,2)
- double precision app_(2,2),bpp_(2,2),rpp_(2,2)
+ double precision app_(2,2),bpp_(2,2),epp_(2,2),rpp_(2,2)
integer ncont,icont(2,maxcont)
double precision econt(maxcont)
double precision boxshift
*
* as of 7/06/91.
*
-c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+c data epp_ / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
zmedi=zi+0.5*dzi
call to_box(xmedi,ymedi,zmedi)
c write (iout,*) "i",xmedi,ymedi,zmedi
-c write (iout,*) "i",xmedi,ymedi,zmedi
do 4 j=i+2,nct-1
c write (iout,*) "i",i," j",j
if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
xj=c(1,j)+0.5*dxj
yj=c(2,j)+0.5*dyj
zj=c(3,j)+0.5*dzj
-c write (iout,*) "j",xj,yj,zj
call to_box(xj,yj,zj)
- xj=boxshift(xj-xi,boxxsize)
- yj=boxshift(yj-yi,boxysize)
- zj=boxshift(zj-zi,boxzsize)
+c write (iout,*) "j",xj,yj,zj
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
+c write (iout,*) "j",xj,yj,zj
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0/(xj*xj+yj*yj+zj*zj)
rmij=sqrt(rrmij)
econt(ncont)=eesij
endif
ees=ees+eesij
-c write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
+c write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
4 continue
1 continue
if (lprint) then
i2=icont(2,i)
it1=itype(i1)
it2=itype(i2)
- write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+ write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)')
& i,restyp(it1),i1,restyp(it2),i2,econt(i)
enddo
endif
i2=icont(2,i)
it1=itype(i1)
it2=itype(i2)
- write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+ write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)')
& i,restyp(it1),i1,restyp(it2),i2,econt(i)
enddo
endif
ii1=max0(ii1-1,1)
jj1=max0(jj1-1,1)
nbeta=nbeta+1
- if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
+ if(lprint)write(iout,'(a,i7,4i6)')'parallel beta',
& nbeta,ii1,i1,jj1,j1
nbfrag=nbfrag+1
if (lprint) then
- write (iout,'(a,i3,4i4)')'antiparallel beta',
+ write (iout,'(a,i3,4i6)')'antiparallel beta',
& nbeta,ii1-1,i1,jj1,j1-1
nstrand=nstrand+1
if (nstrand.le.9) then
- write(12,'(a18,i1,a9,i3,a2,i3,a1)')
+ write(12,'(a18,i1,a9,i6,a2,i6,a1)')
& "DefPropRes 'strand",nstrand,
& "' 'num = ",ii1-2,"..",i1-1,"'"
else
- write(12,'(a18,i2,a9,i3,a2,i3,a1)')
+ write(12,'(a18,i2,a9,i6,a2,i6,a1)')
& "DefPropRes 'strand",nstrand,
& "' 'num = ",ii1-2,"..",i1-1,"'"
endif
nstrand=nstrand+1
if (nstrand.le.9) then
- write(12,'(a18,i1,a9,i3,a2,i3,a1)')
+ write(12,'(a18,i1,a9,i6,a2,i6,a1)')
& "DefPropRes 'strand",nstrand,
& "' 'num = ",j1-2,"..",jj1-1,"'"
else
- write(12,'(a18,i2,a9,i3,a2,i3,a1)')
+ write(12,'(a18,i2,a9,i6,a2,i6,a1)')
& "DefPropRes 'strand",nstrand,
& "' 'num = ",j1-2,"..",jj1-1,"'"
endif
- write(12,'(a8,4i4)')
+ write(12,'(a8,4i6)')
& "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
endif
endif
gg_lipi=0.0d0
gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
gg_lipi=0.0d0
gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
gg_lipi=0.0d0
gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
gg_lipi=0.0d0
gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
c endif
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
c endif
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
& xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision subchap,sss1,sssgrad1
- double precision boxshift
+ double precision boxshift,rij1
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
c if (icall.eq.0) lprn=.false.
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ if (energy_dec)
+ & write(2,*) "g_listscsc_start_long,g_listscsc_end_long",
+ & g_listscsc_start_long,g_listscsc_end_long
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss1=sscale(1.0d0/rij,r_cut_int)
+ rij1=1.0d0/rij
+c sss1=sscale(1.0d0/rij,r_cut_int)
+ sss1=sscale(rij1,r_cut_int)
if (sss1.eq.0.0d0) cycle
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+ rij1=rij1/sigmaii(itypi,itypj)
+ sss=sscale(rij1,r_cut_respa)
+c sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
if (sss.lt.1.0d0) then
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
C assuming the Gay-Berne potential of interaction.
C
implicit none
+ include 'mpif.h'
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
include "COMMON.SPLITELE"
+ include 'COMMON.TIME1'
logical lprn
double precision evdw
integer itypi,itypj,itypi1,iint,ind,ikont
& sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision boxshift
+ double precision time01
+c time01=MPI_Wtime()
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
c if (icall.eq.0) lprn=.false.
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ if (energy_dec)
+ & write(2,*) "g_listscsc_start_short,g_listscsc_end_short",
+ & g_listscsc_start_short,g_listscsc_end_short
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
- sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
if (sss.gt.0.0d0) then
+ sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
c enddo ! j
c enddo ! iint
enddo ! i
+c time_evdw_short=time_evdw_short+MPI_Wtime()-time01
c write (iout,*) "Number of loop steps in EGB:",ind
cccc energy_dec=.false.
return
c if (icall.eq.0) lprn=.true.
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
c if (icall.eq.0) lprn=.true.
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
c do i=iatel_s,iatel_e
+ if (energy_dec)
+ & write(iout,*) "g_listpp_start,g_listpp_end",
+ & g_listpp_start,g_listpp_end
do ikont=g_listpp_start,g_listpp_end
i=newcontlistppi(ikont)
j=newcontlistppj(ikont)
c & " iatel_e_vdw",iatel_e_vdw
c call flush(iout)
c do i=iatel_s_vdw,iatel_e_vdw
- do ikont=g_listpp_vdw_start,g_listpp_vdw_end
- i=newcontlistpp_vdwi(ikont)
- j=newcontlistpp_vdwj(ikont)
+ if (energy_dec)
+ & write(iout,*) "g_listpp_vdw_start_short,g_listpp_vdw_end_short",
+ & g_listpp_vdw_start_short,g_listpp_vdw_end_short
+ do ikont=g_listpp_vdw_start_short,g_listpp_vdw_end_short
+ i=newcontlistpp_vdwi_short(ikont)
+ j=newcontlistpp_vdwj_short(ikont)
if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
c & ' iatscp_e=',iatscp_e
c do i=iatscp_s,iatscp_e
- do ikont=g_listscp_start,g_listscp_end
- i=newcontlistscpi(ikont)
- j=newcontlistscpj(ikont)
+ if (energy_dec)
+ & write(iout,*)"g_listscp_start_long,g_listscp_end_long",
+ & g_listscp_start_long,g_listscp_end_long
+ do ikont=g_listscp_start_long,g_listscp_end_long
+ i=newcontlistscpi_long(ikont)
+ j=newcontlistscpj_long(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
double precision ggg(3)
double precision sscale,sscagrad
double precision boxshift
+ integer ikont
evdw2=0.0D0
evdw2_14=0.0d0
cd print '(a)','Enter ESCP'
c if (lprint_short)
c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
c & ' iatscp_e=',iatscp_e
- if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
- do i=iatscp_s,iatscp_e
+c if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
+ if (energy_dec)
+ & write(iout,*) "g_listscp_start_short,g_listscp_end_short",
+ & g_listscp_start_short,g_listscp_end_short
+ do ikont=g_listscp_start_short,g_listscp_end_short
+ i=newcontlistscpi_short(ikont)
+ j=newcontlistscpj_short(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
c if (lprint_short)
c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
c & " nscp_gr",nscp_gr(i)
- do iint=1,nscp_gr(i)
-
- do j=iscpstart(i,iint),iscpend(i,iint)
+c do iint=1,nscp_gr(i)
+c
+c do j=iscpstart(i,iint),iscpend(i,iint)
itypj=iabs(itype(j))
c if (lprint_short)
c & write (iout,*) "j",j," itypj",itypj
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
endif
- enddo
+c enddo
- enddo ! iint
+c enddo ! iint
enddo ! i
do i=1,nct
do j=1,3
& eliptran,Eafmforce,Etube,
& esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
integer n_corr,n_corr1
+ double precision time01
#ifdef MPI
c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
c & " nfgtasks",nfgtasks
endif
c write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate
if (mod(itime_mat,imatupdate).eq.0) then
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
call make_SCp_inter_list
c write (iout,*) "Finished make_SCp_inter_list"
c call flush(iout)
call make_pp_inter_list
c write (iout,*) "Finished make_pp_inter_list"
c call flush(iout)
- call make_pp_vdw_inter_list
+c call make_pp_vdw_inter_list
c write (iout,*) "Finished make_pp_vdw_inter_list"
c call flush(iout)
+#ifdef TIMING_ENE
+ time_list=time_list+MPI_Wtime()-time01
+#endif
endif
c print *,'Processor',myrank,' calling etotal ipot=',ipot
c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
C Compute the side-chain and electrostatic interaction energy
C
C print *,ipot
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
goto (101,102,103,104,105,106) ipot
C Lennard-Jones potential.
101 call elj(evdw)
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+#ifdef TIMING_ENE
+ time_evdw=time_evdw+MPI_Wtime()-time01
+#endif
#ifdef DFA
C BARTEK for dfa test!
c print *,"Processors",MyRank," wdfa",wdfa_dist
#ifdef TIMING
time_vec=time_vec+MPI_Wtime()-time01
#endif
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
C Introduction of shielding effect first for each peptide group
C the shielding factor is set this factor is describing how each
C peptide group is shielded by side-chains
c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
c & eello_turn4)
endif
+#ifdef TIMING_ENE
+ time_eelec=time_eelec+MPI_Wtime()-time01
+#endif
c#ifdef TIMING
c time_enecalc=time_enecalc+MPI_Wtime()-time00
c#endif
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
if (ipot.lt.6) then
if(wscp.gt.0d0) then
call escp(evdw2,evdw2_14)
c write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
+#ifdef TIMING_ENE
+ time_escp=time_escp+MPI_Wtime()-time01
+#endif
c
c Calculate the bond-stretching energy
c
gradbufc(k,i)=0.0d0
enddo
enddo
-#ifdef DEBUG
- write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
- write (iout,*) (i," jgrad_start",jgrad_start(i),
- & " jgrad_end ",jgrad_end(i),
- & i=igrad_start,igrad_end)
-#endif
+c#ifdef DEBUG
+c write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
+c write (iout,*) (i," jgrad_start",jgrad_start(i),
+c & " jgrad_end ",jgrad_end(i),
+c & i=igrad_start,igrad_end)
+c#endif
c
c Obsolete and inefficient code; we can make the effort O(n) and, therefore,
c do not parallelize this part.
double precision sscale,sscagrad,sscagradlip,sscalelip
double precision gg_lipi(3),gg_lipj(3)
double precision boxshift
+ external boxshift
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
gg_lipi=0.0d0
c write(iout,*) "PRZED ZWYKLE", evdwij
call dyn_ssbond_ene(i,j,evdwij)
c write(iout,*) "PO ZWYKLE", evdwij
+c call flush(iout)
evdw=evdw+evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
-c return
+ return
endif
sigder=-sig*sigsq
c---------------------------------------------------------------
write(iout,*) 'b2=',(b2(k,i-2),k=1,2)
#endif
enddo
- mu=0.0d0
+ mu(:,:nres)=0.0d0
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
C side-chain vectors.
C
implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
+ include 'COMMON.TIME1'
double precision ggg(3)
integer i,iint,j,k,iteli,itypj,subchap,ikont
double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
double precision evdw2,evdw2_14,evdwij
double precision sscale,sscagrad
double precision boxshift
+ external boxshift,to_box
+c#ifdef TIMING_ENE
+c double precision time01
+c#endif
evdw2=0.0D0
evdw2_14=0.0d0
c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
c do i=iatscp_s,iatscp_e
do ikont=g_listscp_start,g_listscp_end
+c#ifdef TIMING_ENE
+c time01=MPI_Wtime()
+c#endif
i=newcontlistscpi(ikont)
j=newcontlistscpj(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
+!DIR$ INLINE
call to_box(xi,yi,zi)
c do iint=1,nscp_gr(i)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
+!DIR$ INLINE
call to_box(xj,yj,zj)
+c#ifdef TIMING_ENE
+c time_escpsetup=time_escpsetup+MPI_Wtime()-time01
+c time01=MPI_Wtime()
+c#endif
+!DIR$ INLINE
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
c print *,xj,yj,zj,'polozenie j'
+c#ifdef TIMING_ENE
+c time_escpsetup=time_escpsetup+MPI_Wtime()-time01
+c time01=MPI_Wtime()
+c#endif
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
c print *,rrij
sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
+c#ifdef TIMING_ENE
+c time_escpcalc=time_escpcalc+MPI_Wtime()-time01
+c#endif
c endif !endif for sscale cutoff
c enddo ! j
if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
& iabs(itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
- ehpb=ehpb+2*eij
+c ehpb=ehpb+2*eij
+ ehpb=ehpb+eij
endif
cd write (iout,*) "eij",eij
cd & ' waga=',waga,' fac=',fac
#ifdef MPI
include "mpif.h"
double precision weights_(n_ene)
- double precision time00,time_Bcast,time_BcastW
+ double precision time00
integer ierror,ierr
#endif
include 'COMMON.SETUP'
include 'COMMON.QRESTR'
include 'COMMON.MD'
include 'COMMON.CONTROL'
+ include 'COMMON.TIME1'
double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
& eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
& escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
& esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
integer i,n_corr,n_corr1
c write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot
+#ifdef TIMING_ENE
+ double precision time01
+#endif
if (modecalc.eq.12.or.modecalc.eq.14) then
#ifdef MPI
c if (fg_rank.eq.0) call int_from_cart1(.false.)
call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
endif
if (mod(itime_mat,imatupdate).eq.0) then
- call make_SCp_inter_list
- call make_SCSC_inter_list
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
+ call make_SCp_inter_list_RESPA
+ call make_SCSC_inter_list_RESPA
call make_pp_inter_list
- call make_pp_vdw_inter_list
+ call make_pp_vdw_inter_list_RESPA
+#ifdef TIMING_ENE
+ time_list=time_list+MPI_Wtime()-time01
+#endif
endif
#endif
cd print *,'nnt=',nnt,' nct=',nct
C
C Compute the side-chain and electrostatic interaction energy
C
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
goto (101,102,103,104,105,106) ipot
C Lennard-Jones potential.
101 call elj_long(evdw)
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+#ifdef TIMING_ENE
+ time_evdw_long=time_evdw_long+MPI_Wtime()-time01
+#endif
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
call vec_and_deriv
+#ifdef TIMING
+ time_vec=time_vec+MPI_Wtime()-time01
+#endif
c write (iout,*) "etotal_long: shield_mode",shield_mode
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
if (shield_mode.eq.1) then
call set_shield_fac
else if (shield_mode.eq.2) then
call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
& eello_turn4)
endif
+#ifdef TIMING_ENE
+ time_eelec_long=time_eelec_long+MPI_Wtime()-time01
+#endif
C
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
if (ipot.lt.6) then
if(wscp.gt.0d0) then
call escp_long(evdw2,evdw2_14)
else
call escp_soft_sphere(evdw2,evdw2_14)
endif
+#ifdef TIMING_ENE
+ time_escp_long=time_escp_long+MPI_Wtime()-time01
+#endif
#ifdef FOURBODY
C
C 12/1/95 Multi-body terms
include 'COMMON.CONTROL'
include 'COMMON.SAXS'
include 'COMMON.TORCNSTR'
+ include 'COMMON.TIME1'
double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
& eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
& escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
& eliptran,Eafmforce,Etube,
& esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
integer i,n_corr,n_corr1
+#ifdef TIMING_ENE
+ double precision time01
+#endif
c write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
c call flush(iout)
if (modecalc.eq.12.or.modecalc.eq.14) then
wsaxs=weights(26)
endif
c write (iout,*),"Processor",myrank," BROADCAST weights"
- call MPI_Bcast(c(1,1),maxres6,MPI_DOUBLE_PRECISION,
+ call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,
& king,FG_COMM,IERR)
c write (iout,*) "Processor",myrank," BROADCAST c"
- call MPI_Bcast(dc(1,1),maxres6,MPI_DOUBLE_PRECISION,
+ call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION,
& king,FG_COMM,IERR)
c write (iout,*) "Processor",myrank," BROADCAST dc"
- call MPI_Bcast(dc_norm(1,1),maxres6,MPI_DOUBLE_PRECISION,
+ call MPI_Bcast(dc_norm(1,1),6*nres,MPI_DOUBLE_PRECISION,
& king,FG_COMM,IERR)
c write (iout,*) "Processor",myrank," BROADCAST dc_norm"
call MPI_Bcast(theta(1),nres,MPI_DOUBLE_PRECISION,
C
C Compute the side-chain and electrostatic interaction energy
C
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
goto (101,102,103,104,105,106) ipot
C Lennard-Jones potential.
101 call elj_short(evdw)
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+#ifdef TIMING_ENE
+ time_evdw_short=time_evdw_short+MPI_Wtime()-time01
+#endif
c
c Calculate the short-range part of Evdwpp
c
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
call evdwpp_short(evdw1)
+#ifdef TIMING_ENE
+ time_eelec_short=time_eelec_short+MPI_Wtime()-time01
+#endif
c
c Calculate the short-range part of ESCp
c
+#ifdef TIMING_ENE
+ time01=MPI_Wtime()
+#endif
if (ipot.lt.6) then
call escp_short(evdw2,evdw2_14)
endif
+#ifdef TIMING_ENE
+ time_escp_short=time_escp_short+MPI_Wtime()-time01
+#endif
c
c Calculate the bond-stretching energy
c
C
C Calculate the SC local energy.
C
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
call vec_and_deriv
+#ifdef TIMING
+ time_vec=time_vec+MPI_Wtime()-time01
+#endif
call esc(escloc)
C
C Calculate the virtual-bond torsional energy.
logical overlap,back,fail
integer nstart
integer i,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
- double precision gen_theta,gen_phi,dist
+ double precision gen_theta,gen_phi,dist,ran_number
cd print *,' CG Processor',me,' maxgen=',maxgen
maxsi=100
cd write (iout,*) 'Gen_Rand_conf: nstart=',nstart
it1=iabs(itype(i-1))
it2=iabs(itype(i-2))
it=iabs(itype(i))
-c print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
-c & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
+ if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
+ vbld(i)=ran_number(3.8d0,10.0d0)
+ vbld_inv(i)=1.0d0/vbld(i)
+ endif
+c write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
+c & ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
phi(i+1)=gen_phi(i+1,it1,it)
if (back) then
phi(i)=gen_phi(i+1,it2,it1)
-c print *,'phi(',i,')=',phi(i)
+c write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
if (it2.ne.10 .and. it2.ne.ntyp1) then
nsi=0
endif
call locate_next_res(i-1)
endif
- theta(i)=gen_theta(it1,phi(i),phi(i+1))
+ if (it1.ne.ntyp1) then
+ theta(i)=gen_theta(it1,phi(i),phi(i+1))
+ else
+ theta(i)=ran_number(1.326d0,2.548d0)
+ endif
+c write (iout,*) "i",i," theta",theta(i)
if (it1.ne.10 .and. it1.ne.ntyp1) then
nsi=0
fail=.true.
call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
nsi=nsi+1
enddo
+c write (iout,*) "After forward SC generation:",nsi,maxsi
if (nsi.gt.maxsi) return1
endif
call locate_next_res(i)
if (overlap(i-1)) then
+c write (iout,*) "overlap",i-1
if (nit.lt.maxnit) then
back=.true.
nit=nit+1
endif
endif
else
+c write (iout,*) "No overlap",i-1
back=.false.
nit=0
i=i+1
double precision function gen_phi(i,it1,it2)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
include "COMMON.TORCNSTR"
include 'COMMON.GEO'
include 'COMMON.BOUNDS'
- if (raw_psipred .or. ndih_constr.eq.0) then
+ include 'COMMON.INTERACT'
+ double precision sumprob(3)
+ double precision pinorm
+ external pinorm
+ if (ndih_constr.eq.0) then
gen_phi=ran_number(-pi,pi)
+ else if (raw_psipred) then
+ if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
+ & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
+ ii=iconstr_dih(i)
+ sumprob(1)=vpsipred(2,ii)
+ sumprob(2)=sumprob(1)+vpsipred(3,ii)
+ sumprob(3)=sumprob(2)+vpsipred(1,ii)
+ aux=ran_number(0.0d0,sumprob(3))
+#ifdef DEBUG
+ write(iout,*)"gen_phi: residue i",i," ii",ii," vpsipred",
+ & vpsipred(2,ii),vpsipred(3,ii),vpsipred(1,ii)," sumprob",
+ & sumprob(1),sumprob(2),sumprob(3)
+ write (iout,*) "aux",aux
+#endif
+ if (aux.le.sumprob(1)) then
+#ifdef DEBUG
+ write (iout,*) "hel:",
+ & (phibound(1,i)-sdihed(1,ndih_constr))*rad2deg,
+ & (phibound(1,i)+sdihed(1,ndih_constr))*rad2deg
+#endif
+ gen_phi=ran_number(phibound(1,i)-sdihed(1,ndih_constr)
+ & ,phibound(1,i)+sdihed(1,ndih_constr))
+ else if (aux.le.sumprob(2)) then
+#ifdef DEBUG
+ write (iout,*) "ext:",
+ & (phibound(2,i)-sdihed(2,ndih_constr))*rad2deg,
+ & (phibound(2,i)+sdihed(2,ndih_constr))*rad2deg
+#endif
+ gen_phi=pinorm(ran_number(phibound(2,i)-sdihed(2,ndih_constr)
+ & ,phibound(2,i)+sdihed(2,ndih_constr)))
+ else
+#ifdef DEBUG
+ write (iout,*) "coil:",-180.0,180.0
+#endif
+ gen_phi=ran_number(-pi,pi)
+ endif
+ else
+ gen_phi=ran_number(-pi,pi)
+ endif
else
C 8/13/98 Generate phi using pre-defined boundaries
gen_phi=ran_number(phibound(1,i),phibound(2,i))
fail=.false.
if (the.eq.0.0D0 .or. the.eq.pi) then
#ifdef MPI
- write (*,'(a,i4,a,i3,a,1pe14.5)')
+ write (iout,'(a,i4,a,i3,a,1pe14.5)')
& 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
#else
cd write (iout,'(a,i3,a,1pe14.5)')
include 'COMMON.IOUNITS'
logical had_overlaps,fail,scfail
integer ioverlap(maxres),ioverlap_last
+ integer maxit_corr /5000/
had_overlaps=.false.
- call overlap_sc_list(ioverlap,ioverlap_last)
+ call overlap_sc_list(ioverlap,ioverlap_last,.false.)
if (ioverlap_last.gt.0) then
write (iout,*) '#OVERLAPing residues ',ioverlap_last
- write (iout,'(18i5)') (ioverlap(k),k=1,ioverlap_last)
+ write (iout,'(15i6)') (ioverlap(k),k=1,ioverlap_last)
had_overlaps=.true.
endif
maxsi=1000
- do k=1,1000
+ do k=1,maxit_corr
if (ioverlap_last.eq.0) exit
do ires=1,ioverlap_last
i=ioverlap(ires)
iti=iabs(itype(i))
- if (iti.ne.10) then
+ if (iti.ne.10 .and. iti.lt.ntyp1) then
nsi=0
fail=.true.
do while (fail.and.nsi.le.maxsi)
if(fail) goto 999
endif
enddo
-
c write (iout,*) "before chaincuild overlap_sc_list: dc0",dc(:,0)
c call chainbuild_extconf
c write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
- call overlap_sc_list(ioverlap,ioverlap_last)
- write (iout,*) 'Overlaping residues ',ioverlap_last,
- & (ioverlap(j),j=1,ioverlap_last)
+ call overlap_sc_list(ioverlap,ioverlap_last,.false.)
+ write (iout,*)'#Overlaping residues @iter',k,":",ioverlap_last
+ write (iout,*)'Residue list:',(ioverlap(j),j=1,ioverlap_last)
enddo
- if (k.le.1000.and.ioverlap_last.eq.0) then
+ if (k.le.maxit_corr.and.ioverlap_last.eq.0) then
scfail=.false.
if (had_overlaps) then
write (iout,*) '#OVERLAPing all corrected after ',k,
else
scfail=.true.
write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
- write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
+ write (iout,'(15i6)') (ioverlap(j),j=1,ioverlap_last)
endif
return
999 continue
- write (iout,'(a30,i5,a12,i4)')
+ write (iout,'(a30,i5,a12,i6)')
& '#OVERLAP FAIL in gen_side after',maxsi,
& 'iter for RES',i
scfail=.true.
return
end
- subroutine overlap_sc_list(ioverlap,ioverlap_last)
- implicit real*8 (a-h,o-z)
+ subroutine overlap_sc_list(ioverlap,ioverlap_last,lprn)
+ implicit none
include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.SETUP"
+ integer ierror
+ integer ioverlap_last_tab(0:max_fg_procs-1),
+ & ioverlap_all(maxres*max_fg_procs),displs(0:max_fg_procs-1)
+#endif
include 'COMMON.GEO'
include 'COMMON.LOCAL'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
include 'COMMON.VAR'
include 'COMMON.CALC'
- logical fail
+ integer ii,itypi,itypj,itypi1,ind,ikont
+ logical fail,lprn
integer ioverlap(maxres),ioverlap_last
- data redfac /0.5D0/
+ double precision redfac /0.5D0/
+ double precision rrij,rij_shift,sig0ij,xi,yi,zi,rcomp,sig
+ double precision dist
- write (iout,*) "overlap_sc_list"
+#ifdef MPI
+ if (nfgtasks.gt.1) then
+ if (fg_rank.eq.0)
+ & call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,king,FG_COMM,
+ & IERROR)
+ call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,king,
+ & FG_COMM,IERROR)
+ call MPI_Bcast(dc_norm(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
+ & king,FG_COMM,IERROR)
+ endif
+#endif
+c write (iout,*) "overlap_sc_list"
c write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
- write(iout,*) "nnt",nnt," nct",nct
+c write(iout,*) "nnt",nnt," nct",nct
ioverlap_last=0
C Check for SC-SC overlaps and mark residues
c print *,'>>overlap_sc nnt=',nnt,' nct=',nct
ind=0
+#ifdef DEBUG
+ write (iout,*) "FG proecssor",fg_rank," g_listscsc_start",
+ & g_listscsc_start," g_listscsc_end",g_listscsc_end
+ write (*,*) "FG proecssor",fg_rank," g_listscsc_start",
+ & g_listscsc_start," g_listscsc_end",g_listscsc_end
+#endif
c do i=iatsc_s,iatsc_e
- do i=nnt,nct
+ do ikont=g_listscsc_start,g_listscsc_end
+c write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
+c & ioverlap_last
+ i=newcontlisti(ikont)
+ j=newcontlistj(ikont)
+c do i=nnt,nct
+c write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
+c & ioverlap_last,"ikont i j",ikont,i,j
itypi=iabs(itype(i))
itypi1=iabs(itype(i+1))
if (itypi.eq.ntyp1) cycle
c
c do iint=1,nint_gr(i)
c do j=istart(i,iint),iend(i,iint)
- do j=i+1,nct
+c do j=i+1,nct
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c if ( rij_shift.le.0.0D0 ) then
if ( rij_shift/sig0ij.le.0.1D0 ) then
c write (iout,*) "overlap",i,j
- write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
- & 'overlap SC-SC: i=',i,' j=',j,
- & ' dist=',dist(nres+i,nres+j),' rcomp=',
- & rcomp,1.0/rij,rij_shift
+ if (lprn) then
+ write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
+ & 'overlap SC-SC: i=',i,' j=',j,
+ & ' dist=',dist(nres+i,nres+j),' rcomp=',
+ & rcomp,1.0/rij,rij_shift
+ write (*,'(a,i2,a,i5,a,i5,a,f10.5,a,3f10.5)')
+ & 'FG processor',fg_rank,' overlap SC-SC: i=',i,' j=',j,
+ & ' dist=',dist(nres+i,nres+j),' rcomp=',
+ & rcomp,1.0/rij,rij_shift
+ endif
ioverlap_last=ioverlap_last+1
ioverlap(ioverlap_last)=i
do k=1,ioverlap_last-1
do k=1,ioverlap_last-1
if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
enddo
+c write(*,*) "FG processor",fg_rank,i,j," ioverlap_last",
+c & ioverlap_last," ioverlap",(ioverlap(k),k=1,ioverlap_last)
endif
- enddo
+c enddo
c enddo
enddo
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "FG Processor",fg_rank," ioverlap_last",
+ & ioverlap_last," ioverlap",(ioverlap(i),i=1,ioverlap_last)
+ write (*,*) "FG Processor",fg_rank," ioverlap_last",ioverlap_last,
+ & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
+ call flush(iout)
+#endif
+ if (nfgtasks.eq.1) return
+#ifdef DEBUG
+ write (iout,*) "Before MPI_Gather"
+ call flush(iout)
+#endif
+ call MPI_Gather(ioverlap_last,1,MPI_INTEGER,ioverlap_last_tab,
+ & 1,MPI_INTEGER,king,FG_COMM,IERROR)
+#ifdef DEBUG
+ write (iout,*) "After MPI_Gather"
+ call flush(iout)
+#endif
+#ifdef DEBUG
+ if (myrank.eq.king)
+ & write (iout,*) "FG Processor",fg_rank,"ioverlap_last_tab",
+ & (ioverlap_last_tab(i),i=0,nfgtasks-1)
+ call flush(iout)
+#endif
+ displs(0)=0
+ do i=1,nfgtasks-1
+ displs(i)=displs(i-1)+ioverlap_last_tab(i-1)
+ enddo
+ call MPI_Gatherv(ioverlap,ioverlap_last,MPI_INTEGER,
+ & ioverlap_all,ioverlap_last_tab,displs,MPI_INTEGER,king,
+ & FG_COMM,IERROR)
+#ifdef DEBUG
+ write (iout,*) "After Gatherv"
+ call flush(iout)
+#endif
+ if (fg_rank.gt.0) return
+ ioverlap_last=0
+ do i=0,nfgtasks-1
+ ioverlap_last=ioverlap_last+ioverlap_last_tab(i)
+ enddo
+#ifdef DEBUG
+ write (iout,*) "ioverlap_last",ioverlap_last," ioverlap_last",
+ & (ioverlap_all(i),i=1,ioverlap_last)
+ call flush(iout)
+#endif
+ ii=0
+ do i=1,ioverlap_last
+ ioverlap(ii+1)=ioverlap_all(i)
+ do j=ii,1,-1
+ if (ioverlap(ii+1).eq.ioverlap(j)) goto 11
+ enddo
+ ii=ii+1
+ 11 continue
+ enddo
+ ioverlap_last=ii
+#ifdef DEBUG
+ write (iout,*) "After summing: ioverlap_last",ioverlap_last,
+ & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
+ call flush(iout)
+#endif
+#endif
return
end
--- /dev/null
+ subroutine gen_rand_conf_mchain(nstart0,*)
+C Generate random conformation or chain cut and regrowth.
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MCM'
+ include 'COMMON.GEO'
+ include 'COMMON.CONTROL'
+ logical overlap_mchain,back,fail
+ integer nstart,nstart0
+ integer i,ii,iii,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
+ integer igr,iequi,ichain,nnres
+ double precision aux
+ double precision gen_theta,gen_phi,dist,ran_number,scalar
+c write (iout,*) 'gen_rand_conf_mchain: maxgen=',maxgen
+ nstart=nstart0
+ maxsi=100
+c write (iout,*) 'Gen_Rand_conf_mchain: nstart=',nstart,
+c & " nchain_group",nchain_group
+
+ DO IGR=1,NCHAIN_GROUP
+
+ DO IEQUI=1,NEQUIV(IGR)
+
+ ichain=iequiv(iequi,igr)
+
+ i=chain_border1(1,ichain)+nstart-1
+ if (nstart.eq.1) then
+ do j=1,3
+ c(j,i)=ran_number(-15.0d0,15.0d0)
+ dc(j,i-1)=c(j,i)
+ enddo
+ endif
+ if (nstart.le.2) then
+
+ do j=1,3
+ dc_norm(j,i)=ran_number(-1.0d0,1.0d0)
+ enddo
+ aux=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
+ do j=1,3
+ dc_norm(j,i)=dc_norm(j,i)/aux
+ enddo
+ if (itype(i).eq.ntyp1) then
+ do j=1,3
+ dc(j,i)=1.9d0*dc_norm(j,i)
+ enddo
+ else
+ do j=1,3
+ dc(j,i)=3.8d0*dc_norm(j,i)
+ enddo
+ endif
+ do j=1,3
+ c(j,i+1)=c(j,i)+dc(j,i)
+ enddo
+ endif
+c if (nstart.lt.5) then
+ if (nstart.le.2) then
+ it1=iabs(itype(i+2))
+ phi(i+3)=gen_phi(i+3,iabs(itype(i+1)),iabs(itype(i+2)))
+c write(iout,*)'phi(4)=',rad2deg*phi(4)
+ theta(i+2)=gen_theta(iabs(itype(i+2)),pi,phi(i+3))
+ if (it1.ne.10) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail)
+ nsi=nsi+1
+ enddo
+ if (nsi.gt.maxsi) then
+ write (iout,'(a,i7,a,i7,a,i7,a,i7)')
+ & 'Problem in SC rotamer generation, residue',ii,
+ & ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+ return1
+ endif
+ endif ! it1.ne.10
+ call orig_frame_chain(i)
+ else
+c call orig_frame_chain(i-1)
+c write (iout,*) "calling refsys",i
+ call refsys(i,i-1,i-2,prod(1,1,i-1),
+ & prod(1,2,i-1),prod(1,3,i-1),fail)
+c write (iout,*) "after refsys",i
+#ifdef DEBUG
+ write (iout,*) "dc_norm(:",i-1,") and prod"
+ do j=1,3
+ write (iout,*) j,dc_norm(j,i-1),(prod(j,k,i-1),k=1,3)
+ enddo
+#endif
+ endif
+
+ ENDDO
+
+ nstart=nstart+1
+ ii=nstart
+
+ maxnit=5000
+
+ nit=0
+ niter=0
+ back=.false.
+ nnres=chain_border1(2,iequiv(1,igr))-
+ & chain_border1(1,iequiv(1,igr))+1
+#ifdef DEBUG
+ write (iout,*) "chain group",igr," chains",
+ & (iequiv(j,igr),j=1,nequiv(igr))
+ write (iout,*) "ii",ii," nnres",nnres
+#endif
+ do while (ii.le. nnres .and. niter.lt.maxgen)
+
+ ichain=iequiv(1,igr)
+ i=ii-1+chain_border1(1,ichain)
+#ifdef DEBUG
+ write (iout,*) "ii",ii," nnres",nnres," ichain",ichain," i",i,
+ & "niter",niter," back",back," nstart",nstart
+#endif
+ if (ii.lt.nstart) then
+c if(iprint.gt.1) then
+ write (iout,'(/80(1h*)/2a/80(1h*))')
+ & 'Generation procedure went down to ',
+ & 'chain beginning. Cannot continue...'
+ write (*,'(/80(1h*)/2a/80(1h*))')
+ & 'Generation procedure went down to ',
+ & 'chain beginning. Cannot continue...'
+c endif
+ return1
+ endif
+ it1=iabs(itype(i-1))
+ it2=iabs(itype(i-2))
+ it=iabs(itype(i))
+ if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
+ vbld(i)=ran_number(3.8d0,10.0d0)
+ vbld_inv(i)=1.0d0/vbld(i)
+ endif
+#ifdef DEBUG
+ write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
+ & ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
+#endif
+ phi(i)=gen_phi(i,it2,it1)
+ phi(i+1)=gen_phi(i+1,it1,it)
+#ifdef DEBUG
+ write (iout,*) "phi",i,phi(i)," phi",i+1,phi(i+1)
+#endif
+ do iequi=2,nequiv(igr)
+ iii=ii+chain_border1(1,iequiv(iequi,igr))
+ phi(iii)=phi(i+1)
+ enddo
+#ifdef CHUJ
+ if (back) then
+ phi(i)=gen_phi(i+1,it2,it1)
+#ifdef DEBUG
+ write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
+#endif
+ theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
+ if (theta(i-1).gt.2.68780478d0) theta(i-1)=2.68780478d0
+ do iequi=2,nequiv(igr)
+ iii=ii-1+chain_border1(1,iequiv(iequi,igr))
+ phi(iii)=phi(i+1)
+ theta(iii-1)=theta(i-1)
+ enddo
+ if (it2.ne.10 .and. it2.ne.ntyp1) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
+ nsi=nsi+1
+ enddo
+ do iequi=2,nequiv(igr)
+ iii=ii-3+chain_border1(1,iequiv(iequi,igr))
+ alph(iii)=alph(i-2)
+ omeg(iii)=omeg(i-2)
+ enddo
+#ifdef DEBUG
+ write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
+#endif
+ if (nsi.gt.maxsi) then
+ write (iout,'(a,i7,a,i7,a,i7,a,i7)')
+ & 'Problem in SC rotamer generation, residue',ii,
+ & ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+ return1
+ endif
+
+ endif
+c call locate_next_res(i-1)
+ endif
+#endif
+ if (it1.ne.ntyp1) then
+ theta(i)=gen_theta(it1,phi(i),phi(i+1))
+ if (theta(i).gt.2.68780478d0) theta(i)=2.68780478d0
+ else
+ theta(i)=ran_number(1.326d0,2.548d0)
+ endif
+#ifdef DEBUG
+ write (iout,*) "ii",ii," i",i," it1",it1,
+ & " theta",theta(i)," phi",
+ & phi(i),phi(i+1)
+#endif
+ if (it1.ne.10 .and. it1.ne.ntyp1) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+c theta(i)=gen_theta(it1,phi(i),phi(i+1))
+ call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
+ nsi=nsi+1
+ enddo
+#ifdef DEBUG
+ write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
+#endif
+c write (iout,*) "After forward SC generation:",nsi,maxsi
+ if (nsi.gt.maxsi) then
+ write (iout,'(a,i7,a,i7,a,i7,a,i7)')
+ & 'Problem in SC rotamer generation, residue',ii,
+ & ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+ return1
+ endif
+
+ endif
+
+ do iequi=2,nequiv(igr)
+ iii=ii-1+chain_border1(1,iequiv(iequi,igr))
+ theta(iii)=theta(i)
+ phi(iii)=phi(i)
+ alph(iii-1)=alph(i-1)
+ omeg(iii-1)=omeg(i-1)
+ enddo
+
+ DO IEQUI=1,NEQUIV(IGR)
+
+ ichain=iequiv(iequi,igr)
+ i=ii-1+chain_border1(1,ichain)
+#ifdef CHUJ
+ if (back) call locate_next_res(i-1)
+#endif
+ call locate_next_res(i)
+#ifdef DEBUG
+ write (iout,*) "i",i," ii",ii," ichain",ichain
+ write (iout,*) theta(i)*rad2deg,phi(i)*rad2deg,
+ & alph(i-1)*rad2deg,omeg(i-1)*rad2deg
+ write (iout,*) (c(j,i),j=1,3),(c(j,i+nres-1),j=1,3)
+#endif
+ if (overlap_mchain(i-1,ii-1,ichain,igr)) then
+#ifdef DEBUG
+ write (iout,*) "***********overlap",i-1," nit",nit
+#endif
+ if (nit.lt.maxnit) then
+ back=.true.
+ nit=nit+1
+ exit
+ else
+#ifdef DEBUG
+ write (iout,*) "***********overlap maxnit exceeded",nit
+#endif
+ nit=0
+ if (ii.gt.3) then
+ back=.true.
+ ii=ii-1
+c write (iout,*) "ii",ii
+ exit
+ else
+ write (iout,'(a,i7,a,i7,a,i7,a,i7)')
+ & 'Cannot generate non-overlaping conformation, residue',ii,
+ & ' chain',ichain,' chain group',igr,'. Increase MAXNIT.'
+ return1
+ endif
+ endif
+ else
+c write (iout,*) "No overlap",i-1
+ back=.false.
+c nit=0
+c i=i+1
+ endif
+
+ ENDDO
+
+ if (.not.back) then
+#ifdef DEBUG
+ write (iout,*) "++++++++++Successful generation",igr,ichain,ii
+#endif
+ ii=ii+1
+ nit=0
+ endif
+ back=.false.
+c write (iout,*) "ii",ii
+ niter=niter+1
+
+ enddo
+ if (niter.ge.maxgen) then
+ write (iout,'(a,2i7,a,i7,a,i7,a,i7)')
+ & 'Too many trials in conformation generation',niter,maxgen,
+ & ' chain group',igr,' chain',ichain,' residue',ii
+ return1
+ endif
+
+ ENDDO
+
+ do i=2,nres
+ if (itype(i).eq.ntyp1) then
+ do j=1,3
+ dc(j,i-1)=c(j,i)-c(j,i-1)
+ enddo
+ endif
+ enddo
+
+ return
+ end
+c-------------------------------------------------------------------------
+ logical function overlap_mchain(i,ii,ichain,igr)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ double precision redfac /0.5D0/,redfacp /0.8d0/,redfacscp /0.8d0/
+ integer i,ii,ichain,igr,j,jj,jchain,jgr,k,iti,itj,iteli,itelj
+ double precision rcomp
+ double precision dist
+ logical lprn /.false./
+ overlap_mchain=.false.
+ iti=iabs(itype(i))
+ if (iti.gt.ntyp) return
+C Check for SC-SC overlaps.
+cd print *,'nnt=',nnt,' nct=',nct
+#ifdef DEBUG
+ write (iout,*) "overlap_mchain i",i," ii",ii," ichain",ichain,
+ & " igr",igr," iti",iti
+#endif
+ do j=nnt,nct
+ if (itype(j).eq.ntyp1) cycle
+ jchain=ireschain(j)
+c write (iout,*) "overlap_mchain j",j," jj",jj," jchain",jchain,
+c & " itype",itype(j)
+ jgr=mapchain(jchain)
+ jj=j-chain_border1(1,jchain)+1
+ itj=iabs(itype(j))
+c write (iout,*) "jgr",jgr
+ if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-1
+ & .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+ if (j.lt.i-1 .or. ipot.ne.4) then
+ rcomp=sigmaii(iti,itj)
+ else
+ rcomp=sigma(iti,itj)
+ endif
+cd print *,'j=',j
+ if (dist(nres+i,nres+j).lt.redfac*rcomp) then
+ overlap_mchain=.true.
+ if (lprn) write(iout,*)'overlap_mchain, SC-SC: i=',i,' j=',j,
+ & ' ichain',ichain,' jchain',jchain,
+ & ' dist=',dist(nres+i,nres+j),' rcomp=',
+ & rcomp*redfac
+ return
+ endif
+ enddo
+#ifdef CHUJ
+C Check for overlaps between the added peptide group and the preceding
+C SCs.
+ iteli=itel(i)
+ if (iteli.gt.0) then
+ do j=1,3
+ c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+ enddo
+ do j=nnt,nct
+ itj=iabs(itype(j))
+ if (itj.eq.ntyp1) cycle
+ jchain=ireschain(j)
+ jgr=mapchain(jchain)
+ jj=j-chain_border1(1,jchain)+1
+ if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
+ & .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+cd print *,'overlap_mchain, p-Sc: i=',i,' j=',j,
+cd & ' dist=',dist(nres+j,maxres2+1)
+ if (dist(j,maxres2+1).lt.4.0D0*redfacscp) then
+ if (lprn) write (iout,*) 'overlap_mchain, p-SC: i=',i,' j=',j,
+ & ' ichain',ichain,' jchain',jchain,
+ & ' dist=',dist(nres+j,maxres2+1),' rcomp=',
+ & 4.0d0*redfac
+ overlap_mchain=.true.
+ return
+ endif
+ enddo
+ endif
+C Check for p-p overlaps
+#endif
+ iteli=itel(i)
+ if (iteli.eq.0) return
+ do j=1,3
+ c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+ enddo
+c do j=nnt,i-2
+ do j=nnt,nct
+ itelj=itel(j)
+ if (itelj.eq.0) cycle
+ do k=1,3
+ c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
+ enddo
+ jchain=ireschain(j)
+ jgr=mapchain(jchain)
+c if(iteli.ne.0.and.itelj.ne.0)then
+c write (iout,*) i,j,dist(maxres2+1,maxres2+2),rpp(iteli,itelj)
+ jj=j-chain_border1(1,jchain)+1
+c write (iout,*) "jgr",jgr
+ if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
+ & .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+#ifdef DEBUG
+ write (iout,*)'overlap_mchain, p-p: i=',i,' j=',j,' igr',igr,
+ & ' jgr',jgr,' ichain',ichain,' jchain',jchain,
+ & ' dist=',dist(maxres2+1,maxres2+2)
+#endif
+ if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfacp) then
+ if (lprn) write (iout,*) 'overlap_mchain, p-p: i=',i,' j=',j,
+ & ' ichain',ichain,' jchain',jchain,
+ & ' dist=',dist(maxres2+1,maxres2+2),' rcomp=',
+ & rpp(iteli,itelj)*redfacp
+ overlap_mchain=.true.
+ return
+ endif
+c endif
+ enddo
+ return
+ end
+
#endif
character*50 tytul
integer iunit
- character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/
+ character*1 chainid(52) /'A','B','C','D','E','F','G','H','I','J',
+ & 'K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
+ & 'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p',
+ & 'q','r','s','t','u','v','w','x','y','z'/
integer ica(maxres)
integer i,j,k,iti,itj,itk,itl,iatom,ichain,ires
double precision etot
do i=1,nss
if (dyn_ss) then
write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)')
- & 'SSBOND',i,'CYS',idssb(i)-nnt+1,
- & 'CYS',jdssb(i)-nnt+1
+ & 'SSBOND',i,'CYS',iss(idssb(i))-nnt+1,
+ & 'CYS',iss(jdssb(i))-nnt+1
else
write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)')
& 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,
iti=itype(i)
if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then
ichain=ichain+1
+ if (ichain.gt.52) ichain=1
ires=0
write (iunit,'(a)') 'TER'
- else
+ else if (iti.ne.ntyp1) then
ires=ires+1
iatom=iatom+1
ica(i)=iatom
write (iunit,'(a6)') 'ENDMDL'
10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
- 30 FORMAT ('CONECT',8I5)
+ 30 FORMAT ('CONECT',8I7)
return
end
c------------------------------------------------------------------------------
if(itime.eq.0) then
write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
& ",31a12)"
- write (istat,format) "#","",
+ write (istat,format) "#"," ",
& (ename(print_order(i)),i=1,nprint_ene)
endif
write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
include 'COMMON.IOUNITS'
include 'COMMON.TIME1'
integer i,j,kk
+ double precision time00,time01
c
c This subrouting calculates total Cartesian coordinate gradient.
c The subroutine chainbuild_cart and energy MUST be called beforehand.
1 "WSC ","WSCP ","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
! 8 9 10 11 12 13 14
8 "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR ","WTORD ",
-! 15 16 17 18 19 20 21
+! 15 16 17 18 19 20 21
5 "WSTRAIN","WVDWPP","WBOND","SCAL14","WDIHC","WUMB","WSCCOR",
! 22 23 24 25 26 27 28
2 "WLT","WAFM","WTHETCNSR","WTUBE","WSAXS","WHOMO","WDFAD",
integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
& iturn4_end_all,iatel_s_all,
& iatel_e_all,ielstart_all,ielend_all,ntask_cont_from_all,
- & itask_cont_from_all,ntask_cont_to_all,itask_cont_to_all,
- & n_sc_int_tot,my_sc_inds,my_sc_inde,ind_sctint,ind_scint_old
+ & itask_cont_from_all,ntask_cont_to_all,itask_cont_to_all
+ integer*8 n_sc_int_tot,my_sc_inds,my_sc_inde,ind_scint,
+ & ind_scint_old,nele_int_tot,ind_eleint,my_ele_inds,my_ele_inde,
+ & ind_eleint_old,nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw,
+ & ind_eleint_vdw,ind_eleint_vdw_old,nscp_int_tot,my_scp_inds,
+ & my_scp_inde,ind_scpint,ind_scpint_old,nsumgrad,nlen,ngrad_start,
+ & ngrad_end
common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
& iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
& iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
& itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
logical scheck,lprint,flag
- integer i,j,k,ii,jj,iint,npept,nele_int_tot,ind_eleint,ind_scint,
- & my_ele_inds,my_ele_inde,ind_eleint_old,nele_int_tot_vdw,
- & my_ele_inds_vdw,my_ele_inde_vdw,ind_eleint_vdw,ijunk,
- & ind_eleint_vdw_old,nscp_int_tot,my_scp_inds,my_scp_inde,
- & ind_scpint,ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,
- & iaux,ind_typ,ncheck_from,ncheck_to,ichunk
+ integer i,j,k,ii,jj,iint,npept,
+ & ijunk,iaux,ind_typ,ncheck_from,ncheck_to,ichunk
#ifdef MPI
integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
& my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
lprint=energy_dec
if (lprint)
&write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
- n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
- call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
+ n_sc_int_tot=int(nct-nnt+1,8)*int(nct-nnt,8)/2-nss
+ call int_bounds8(n_sc_int_tot,my_sc_inds,my_sc_inde)
if (lprint)
& write (iout,*) 'Processor',fg_rank,' CG group',kolor,
& ' absolute rank',MyRank,
& ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
& ' my_sc_inde',my_sc_inde
- ind_sctint=0
+ ind_scint=0
iatsc_s=0
iatsc_e=0
#endif
if (jj.eq.i+1) then
#ifdef MPI
c write (iout,*) 'jj=i+1'
- call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+ call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
& iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
#else
nint_gr(i)=1
else if (jj.eq.nct) then
#ifdef MPI
c write (iout,*) 'jj=nct'
- call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+ call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
& iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
#else
nint_gr(i)=1
#endif
else
#ifdef MPI
- call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+ call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
& iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
ii=nint_gr(i)+1
- call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+ call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
& iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
#else
nint_gr(i)=2
endif
else
#ifdef MPI
- call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+ call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
& iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
#else
nint_gr(i)=1
& ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
#endif
if (lprint) then
+ write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
write (iout,'(a)') 'Interaction array:'
do i=iatsc_s,iatsc_e
- write (iout,'(i3,2(2x,2i3))')
+ write (iout,'(i7,2(2x,2i7))')
& i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
enddo
endif
#ifdef MPI
C Now partition the electrostatic-interaction array
npept=nct-nnt
- nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
- call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
+ nele_int_tot=int(npept-ispp,8)*int(npept-ispp+1,8)/2
+ call int_bounds8(nele_int_tot,my_ele_inds,my_ele_inde)
if (lprint)
& write (*,*) 'Processor',fg_rank,' CG group',kolor,
& ' absolute rank',MyRank,
ind_eleint_old=0
do i=nnt,nct-3
ijunk=0
- call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
+ call int_partition8(ind_eleint,my_ele_inds,my_ele_inde,i,
& iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
enddo ! i
13 continue
if (iatel_s.eq.0) iatel_s=1
- nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
+ nele_int_tot_vdw=int(npept-2,8)*int(npept-2+1,8)/2
c write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
- call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
+ call int_bounds8(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
c write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
c & " my_ele_inde_vdw",my_ele_inde_vdw
ind_eleint_vdw=0
iatel_e_vdw=0
do i=nnt,nct-3
ijunk=0
- call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
+ call int_partition8(ind_eleint_vdw,my_ele_inds_vdw,
& my_ele_inde_vdw,i,
& iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
& ielend_vdw(i),*15)
& ' absolute rank',MyRank
write (iout,*) 'Electrostatic interaction array:'
do i=iatel_s,iatel_e
- write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
+ write (iout,'(i7,2(2x,2i7))') i,ielstart(i),ielend(i)
enddo
endif ! lprint
c iscp=3
iscp=2
C Partition the SC-p interaction array
#ifdef MPI
- nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
- call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
+ nscp_int_tot=int(npept-iscp+1,8)*int(npept-iscp+1,8)
+ call int_bounds8(nscp_int_tot,my_scp_inds,my_scp_inde)
if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
& ' absolute rank',myrank,
& ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
do i=nnt,nct-1
if (i.lt.nnt+iscp) then
cd write (iout,*) 'i.le.nnt+iscp'
- call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+ call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
& iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
& iscpend(i,1),*14)
else if (i.gt.nct-iscp) then
cd write (iout,*) 'i.gt.nct-iscp'
- call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+ call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
& iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
& iscpend(i,1),*14)
else
- call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+ call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
& iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
& iscpend(i,1),*14)
ii=nscp_gr(i)+1
- call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+ call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
& iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
& iscpend(i,ii),*14)
endif
if (lprint) then
write (iout,'(a)') 'SC-p interaction array:'
do i=iatscp_s,iatscp_e
- write (iout,'(i3,2(2x,2i3))')
+ write (iout,'(i7,2(2x,2i7))')
& i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
enddo
endif ! lprint
endif
c nsumgrad=(nres-nnt)*(nres-nnt+1)/2
c nlen=nres-nnt+1
- nsumgrad=(nres-nnt)*(nres-nnt+1)/2
- nlen=nres-nnt+1
- call int_bounds(nsumgrad,ngrad_start,ngrad_end)
- igrad_start=((2*nlen+1)
- & -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
- jgrad_start(igrad_start)=
- & ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
- & +igrad_start
- jgrad_end(igrad_start)=nres
- igrad_end=((2*nlen+1)
- & -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
- if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
- jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
- & +igrad_end
- do i=igrad_start+1,igrad_end-1
- jgrad_start(i)=i+1
- jgrad_end(i)=nres
- enddo
+c nsumgrad=(nres-nnt)*(nres-nnt+1)/2
+c nlen=nres-nnt+1
+c call int_bounds(nsumgrad,ngrad_start,ngrad_end)
+c igrad_start=((2*nlen+1)
+c & -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
+c jgrad_start(igrad_start)=
+c & ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
+c & +igrad_start
+c jgrad_end(igrad_start)=nres
+c igrad_end=((2*nlen+1)
+c & -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
+c if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
+c jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
+c & +igrad_end
+c do i=igrad_start+1,igrad_end-1
+c jgrad_start(i)=i+1
+c jgrad_end(i)=nres
+c enddo
if (lprint) then
write (*,*) 'Processor:',fg_rank,' CG group',kolor,
& ' absolute rank',myrank,
& ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
& ithetaconstr_end
- write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
- & igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
- & ' ngrad_end',ngrad_end
- do i=igrad_start,igrad_end
- write(*,*) 'Processor:',fg_rank,myrank,i,
- & jgrad_start(i),jgrad_end(i)
- enddo
+c write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
+c & igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
+c & ' ngrad_end',ngrad_end
+c do i=igrad_start,igrad_end
+c write(*,*) 'Processor:',fg_rank,myrank,i,
+c & jgrad_start(i),jgrad_end(i)
+c enddo
endif
if (nfgtasks.gt.1) then
call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
return
end
c---------------------------------------------------------------------------
+ subroutine int_bounds8(total_ints,lower_bound,upper_bound)
+ implicit none
+ include 'DIMENSIONS'
+ include 'mpif.h'
+ include 'COMMON.SETUP'
+ integer*8 total_ints,lower_bound,upper_bound,nint
+ integer*8 int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
+ integer i,nexcess
+ nint=total_ints/nfgtasks
+ do i=1,nfgtasks
+ int4proc(i-1)=nint
+ enddo
+ nexcess=total_ints-nint*nfgtasks
+ do i=1,nexcess
+ int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
+ enddo
+ lower_bound=0
+ do i=0,fg_rank-1
+ lower_bound=lower_bound+int4proc(i)
+ enddo
+ upper_bound=lower_bound+int4proc(fg_rank)
+ lower_bound=lower_bound+1
+ return
+ end
+c---------------------------------------------------------------------------
subroutine int_bounds1(total_ints,lower_bound,upper_bound)
implicit none
include 'DIMENSIONS'
endif
return
end
+c---------------------------------------------------------------------------
+ subroutine int_partition8(int_index,lower_index,upper_index,atom,
+ & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ integer*8 int_index,lower_index,upper_index
+ integer atom,at_start,at_end,
+ & first_atom,last_atom,int_gr,jat_start,jat_end,int_index_old
+ logical lprn
+ lprn=.false.
+ if (lprn) write (iout,*) 'int_index=',int_index
+ int_index_old=int_index
+ int_index=int_index+last_atom-first_atom+1
+ if (lprn)
+ & write (iout,*) 'int_index=',int_index,
+ & ' int_index_old',int_index_old,
+ & ' lower_index=',lower_index,
+ & ' upper_index=',upper_index,
+ & ' atom=',atom,' first_atom=',first_atom,
+ & ' last_atom=',last_atom
+ if (int_index.ge.lower_index) then
+ int_gr=int_gr+1
+ if (at_start.eq.0) then
+ at_start=atom
+ jat_start=first_atom-1+lower_index-int_index_old
+ else
+ jat_start=first_atom
+ endif
+ if (lprn) write (iout,*) 'jat_start',jat_start
+ if (int_index.ge.upper_index) then
+ at_end=atom
+ jat_end=first_atom-1+upper_index-int_index_old
+ return1
+ else
+ jat_end=last_atom
+ endif
+ if (lprn) write (iout,*) 'jat_end',jat_end
+ endif
+ return
+ end
#endif
c------------------------------------------------------------------------------
subroutine hpb_partition
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
- integer time00
+ double precision time00
#endif
include 'COMMON.VAR'
include 'COMMON.CHAIN'
#endif
#ifdef FIVEDIAG
call grad_transform
- d_a=0.0d0
+ d_a(:,:2*nres)=0.0d0
if (lprn) then
write (iout,*) "Potential forces backbone"
do i=1,nres
endif
enddo
enddo
- DMorig=DM
- DU1orig=DU1
- DU2orig=DU2
+ DMorig(:2*nres)=DM(:2*nres)
+ DU1orig(:2*nres)=DU1(:2*nres)
+ DU2orig(:2*nres)=DU2(:2*nres)
if (gmatout) then
write (iout,*)"The upper part of the five-diagonal inertia matrix"
endif
c---------------------------------------------------------------------------
subroutine fivediaginv_mult(ndim,forces,d_a_vec)
implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
include 'DIMENSIONS'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
& xsolv(ndim),d_a_vec(6*nres)
integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
- accel=0.0d0
+#ifdef TIMING
+ include 'COMMON.TIME1'
+ double precision time01
+ time01=MPI_Wtime()
+#endif
+ accel(:,:2*nres)=0.0d0
do j=1,3
Compute accelerations in Calpha and SC
do ichain=1,nchain
ind=ind+3
endif
enddo
+#ifdef TIMING
+ time_ginvmult=time_ginvmult+MPI_Wtime()-time01
+#endif
#ifdef DEBUG
write (iout,*) "d_a_vec"
write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
! print *,"START make_SC"
#ifdef DEBUG
write (iout,*) "make_SCSC_inter_list maxint_res",maxint_res
+ write (iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
#endif
r_buff_list=5.0d0
ilist_sc=0
#ifdef MPI
#ifdef DEBUG
write (iout,*) "before MPIREDUCE",ilist_sc
- do i=1,ilist_sc
- write (iout,*) i,contlisti(i),contlistj(i)
- enddo
+c do i=1,ilist_sc
+c write (iout,*) i,contlisti(i),contlistj(i)
+c enddo
#endif
if (nfgtasks.gt.1)then
write (iout,*) i,newcontlisti(i),newcontlistj(i)
enddo
#endif
- call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+ call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+#ifdef DEBUG
+ write (iout,*) "g_listscsc_start",g_listscsc_start,
+ & "g_listscsc_end",g_listscsc_end
return
+#endif
end
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine make_SCSC_inter_list_RESPA
+ implicit none
+ include "DIMENSIONS"
+#ifdef MPI
+ include 'mpif.h'
+ include "COMMON.SETUP"
+#endif
+ include "COMMON.CONTROL"
+ include "COMMON.CHAIN"
+ include "COMMON.INTERACT"
+ include "COMMON.SPLITELE"
+ include "COMMON.IOUNITS"
+ double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp
+ double precision dist_init, dist_temp,r_buff_list
+ integer contlisti_long(maxint_res*maxres),
+ & contlisti_short(maxint_res*maxres),
+ & contlistj_long(maxint_res*maxres),
+ & contlistj_short(maxint_res*maxres)
+! integer :: newcontlisti(200*nres),newcontlistj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+ & ilist_sc_long,g_ilist_sc_long,ilist_sc_short,g_ilist_sc_short
+ integer displ(0:max_fg_procs),i_ilist_sc_long(0:max_fg_procs),
+ & i_ilist_sc_short(0:max_fg_procs),ierr
+ logical lprn /.false./
+ double precision boxshift
+ double precision d_scale,r_respa_buf
+! print *,"START make_SC"
+#ifdef DEBUG
+ write (iout,*) "make_SCSC_inter_list maxint_res",maxint_res
+ write (iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
+#endif
+ r_buff_list=5.0d0
+ r_respa_buf=rlamb
+ ilist_sc_long=0
+ ilist_sc_short=0
+ do i=iatsc_s,iatsc_e
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dist_init=dsqrt(xj*xj+yj*yj+zj*zj)
+! r_buff_list is a read value for a buffer
+ if (dist_init.le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+ d_scale=dist_init/sigmaii(itypi,itypj)
+ if (d_scale.le.r_cut_respa+r_respa_buf) then
+ ilist_sc_short=ilist_sc_short+1
+ contlisti_short(ilist_sc_short)=i
+ contlistj_short(ilist_sc_short)=j
+ endif
+ if (d_scale.gt.r_cut_respa-rlamb-r_respa_buf) then
+ ilist_sc_long=ilist_sc_long+1
+! this can be substituted by cantor and anti-cantor
+ contlisti_long(ilist_sc_long)=i
+ contlistj_long(ilist_sc_long)=j
+ endif
+ endif
+ enddo
+ enddo
+ enddo
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE ilist_sc_long",ilist_sc_long
+c do i=1,ilist_sc_long
+c write (iout,*) i,contlisti_long(i),contlistj_long(i)
+c enddo
+ write (iout,*) "before MPIREDUCE ilist_sc_short",ilist_sc_short
+c do i=1,ilist_sc_short
+c write (iout,*) i,contlisti_short(i),contlistj_short(i)
+c enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_sc_long,g_ilist_sc_long,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+ call MPI_Reduce(ilist_sc_short,g_ilist_sc_short,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c write (iout,*) "SCSC after reduce ierr",ierr
+ if (fg_rank.eq.0.and.(g_ilist_sc_long.gt.maxres*maxint_res .or.
+ & g_ilist_sc_short.gt.maxres*maxint_res)) then
+ if ((me.eq.king.or.out1file).and.energy_dec) then
+ write (iout,*) "Too many SCSC interactions",
+ & g_ilist_sc_long,g_ilist_sc_short,
+ & " only",maxres*maxint_res," allowed."
+ write (iout,*) "Specify a smaller r_cut_int and resubmit"
+ call flush(iout)
+ endif
+ write (*,*) "Processor:",me,": Too many SCSC interactions",
+ & g_ilist_sc_long+g_ilist_sc_short," only",
+ & maxres*maxint_res," allowed."
+ write (*,*) "Specify a smaller r_cut_int and resubmit"
+ call MPI_Abort(MPI_COMM_WORLD,ierr)
+ endif
+c write(iout,*) "before bcast",g_ilist_sc_long
+ call MPI_Gather(ilist_sc_long,1,MPI_INTEGER,
+ & i_ilist_sc_long,1,MPI_INTEGER,king,FG_COMM,IERR)
+c write (iout,*) "SCSC after gather ierr",ierr
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_sc_long(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlisti_long,ilist_sc_long,MPI_INTEGER,
+ & newcontlisti_long,i_ilist_sc_long,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCSC after gatherv ierr",ierr
+ call MPI_Gatherv(contlistj_long,ilist_sc_long,MPI_INTEGER,
+ & newcontlistj_long,i_ilist_sc_long,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_sc_long,1,MPI_INT,king,FG_COMM,IERR)
+c write (iout,*) "SCSC bcast reduce ierr",ierr
+! write(iout,*) "before bcast",g_ilist_sc_long
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlisti_long,g_ilist_sc_long,MPI_INT,king,
+ & FG_COMM,IERR)
+c write (iout,*) "SCSC bcast reduce ierr",ierr
+ call MPI_Bcast(newcontlistj_long,g_ilist_sc_long,MPI_INT,king,
+ & FG_COMM,IERR)
+c write (iout,*) "SCSC after bcast ierr",ierr
+! write(iout,*) "before gather",displ(0),displ(1)
+c write(iout,*) "before bcast",g_ilist_sc_short
+ call MPI_Gather(ilist_sc_short,1,MPI_INTEGER,
+ & i_ilist_sc_short,1,MPI_INTEGER,king,FG_COMM,IERR)
+c write (iout,*) "SCSC after gather ierr",ierr
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_sc_short(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlisti_short,ilist_sc_short,MPI_INTEGER,
+ & newcontlisti_short,i_ilist_sc_short,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCSC after gatherv ierr",ierr
+ call MPI_Gatherv(contlistj_short,ilist_sc_short,MPI_INTEGER,
+ & newcontlistj_short,i_ilist_sc_short,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_sc_short,1,MPI_INT,king,FG_COMM,IERR)
+c write (iout,*) "SCSC bcast reduce ierr",ierr
+! write(iout,*) "before bcast",g_ilist_sc_short
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlisti_short,g_ilist_sc_short,MPI_INT,king,
+ & FG_COMM,IERR)
+c write (iout,*) "SCSC bcast reduce ierr",ierr
+ call MPI_Bcast(newcontlistj_short,g_ilist_sc_short,MPI_INT,king,
+ ^ FG_COMM,IERR)
+c write (iout,*) "SCSC after bcast ierr",ierr
+ else
+#endif
+ g_ilist_sc_long=ilist_sc_long
+
+ do i=1,ilist_sc_long
+ newcontlisti_long(i)=contlisti_long(i)
+ newcontlistj_long(i)=contlistj_long(i)
+ enddo
+
+ g_ilist_sc_short=ilist_sc_short
+
+ do i=1,ilist_sc_short
+ newcontlisti_short(i)=contlisti_short(i)
+ newcontlistj_short(i)=contlistj_short(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+ if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+ & write (iout,'(a30,2i10,a,2i4)')
+ & "Number of long- and short-range SC-SC interactions",
+ & g_ilist_sc_long,g_ilist_sc_short," per residue on average",
+ & g_ilist_sc_long/nres,g_ilist_sc_short/nres
+#ifdef DEBUG
+ write (iout,*)
+ & "make_SCSC_inter_list: g_ilist_sc_long after GATHERV",
+ & g_ilist_sc_long
+ write (iout,*) "List of long-range SCSC interactions"
+ do i=1,g_ilist_sc_long
+ write (iout,*) i,newcontlisti_long(i),newcontlistj_long(i)
+ enddo
+ write (iout,*)
+ & "make_SCSC_inter_list: g_ilist_sc_short after GATHERV",
+ & g_ilist_sc_short
+ write (iout,*) "List of short-range SCSC interactions"
+ do i=1,g_ilist_sc_short
+ write (iout,*) i,newcontlisti_short(i),newcontlistj_short(i)
+ enddo
+#endif
+ call int_bounds(g_ilist_sc_long,g_listscsc_start_long,
+ & g_listscsc_end_long)
+ call int_bounds(g_ilist_sc_short,g_listscsc_start_short,
+ & g_listscsc_end_short)
+#ifdef DEBUG
+ write (iout,*) "g_list_sc_start",g_listscsc_start_long,
+ & "g_list_sc_end",g_listscsc_end_long
+ write (iout,*)"g_list_sc_start_short",g_listscsc_start_short,
+ & "g_list_sc_end_short",g_listscsc_end_short
+#endif
+ return
+ end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine make_SCp_inter_list
implicit none
include "DIMENSIONS"
write (iout,*) "make_SCp_inter_list maxint_res",maxint_res
#endif
r_buff_list=5.0
- ilist_scp=0
- ilist_scp_first=0
+ ilist_scp=0
+ ilist_scp_first=0
do i=iatscp_s,iatscp_e
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
xi=0.5D0*(c(1,i)+c(1,i+1))
! enddo
#endif
call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+#ifdef DEBUG
+ write (iout,*) "g_listscp_start",g_listscp_start,
+ & "g_listscp_end",g_listscp_end
+#endif
+ return
+ end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine make_SCp_inter_list_RESPA
+ implicit none
+ include "DIMENSIONS"
+#ifdef MPI
+ include 'mpif.h'
+ include "COMMON.SETUP"
+#endif
+ include "COMMON.CONTROL"
+ include "COMMON.CHAIN"
+ include "COMMON.INTERACT"
+ include "COMMON.SPLITELE"
+ include "COMMON.IOUNITS"
+ double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp
+ double precision dist_init, dist_temp,r_buff_list
+ integer contlistscpi_long(2*maxint_res*maxres),
+ & contlistscpi_short(2*maxint_res*maxres),
+ & contlistscpj_long(2*maxint_res*maxres),
+ & contlistscpj_short(2*maxint_res*maxres)
+! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
+ integer i,j,iteli,itypj,subchap,xshift,yshift,zshift,iint,
+ & ilist_scp_long,ilist_scp_short,g_ilist_scp_long,g_ilist_scp_short
+ integer displ(0:max_fg_procs),i_ilist_scp_long(0:max_fg_procs),
+ & i_ilist_scp_short(0:max_fg_procs),ierr
+c integer contlistscpi_f(2*maxint_res*maxres),
+c & contlistscpj_f(2*maxint_res*maxres)
+ double precision boxshift
+ double precision d_scale,r_respa_buf
+! print *,"START make_SC"
+#ifdef DEBUG
+ write (iout,*) "make_SCp_inter_list maxint_res",maxint_res
+#endif
+ r_buff_list=5.0
+ r_respa_buf=rlamb
+ ilist_scp_long=0
+ ilist_scp_short=0
+ do i=iatscp_s,iatscp_e
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ xi=0.5D0*(c(1,i)+c(1,i+1))
+ yi=0.5D0*(c(2,i)+c(2,i+1))
+ zi=0.5D0*(c(3,i)+c(3,i+1))
+ call to_box(xi,yi,zi)
+ iteli=itel(i)
+ do iint=1,nscp_gr(i)
+ do j=iscpstart(i,iint),iscpend(i,iint)
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
+! Uncomment following three lines for SC-p interactions
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+! Uncomment following three lines for Ca-p interactions
+! xj=c(1,j)-xi
+! yj=c(2,j)-yi
+! zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dist_init=dsqrt(xj*xj+yj*yj+zj*zj)
+! r_buff_list is a read value for a buffer
+ if (dist_init.le.(r_cut_int+r_buff_list)) then
+ d_scale=dist_init/rscp(itypj,iteli)
+ if (d_scale.le.r_cut_respa+r_respa_buf) then
+! Here the list is created
+ ilist_scp_short=ilist_scp_short+1
+ contlistscpi_short(ilist_scp_short)=i
+ contlistscpj_short(ilist_scp_short)=j
+ endif
+ if (d_scale.gt.r_cut_respa-rlamb-r_respa_buf) then
+! this can be substituted by cantor and anti-cantor
+ ilist_scp_long=ilist_scp_long+1
+ contlistscpi_long(ilist_scp_long)=i
+ contlistscpj_long(ilist_scp_long)=j
+ endif
+ endif
+ enddo
+ enddo
+ enddo
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_scp_long,ilist_scp_short
+ write (iout,*) "Long-range scp interaction list"
+ do i=1,ilist_scp_long
+ write (iout,*) i,contlistscpi_long(i),contlistscpj_long(i)
+ enddo
+ write (iout,*) "Short-range scp interaction list"
+ do i=1,ilist_scp_short
+ write (iout,*) i,contlistscpi_short(i),contlistscpj_short(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_scp_long,g_ilist_scp_long,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+ call MPI_Reduce(ilist_scp_short,g_ilist_scp_short,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c write (iout,*) "SCp after reduce ierr",ierr
+ if (fg_rank.eq.0.and.(g_ilist_scp_long.gt.
+ & 2*maxres*maxint_res .or. g_ilist_scp_short.gt.
+ & 2*maxres*maxint_res)) then
+ if ((me.eq.king.or.out1file).and.energy_dec) then
+ write (iout,*) "Too many SCp interactions",
+ & g_ilist_scp_long+g_ilist_scp_short," only",
+ & 2*maxres*maxint_res," allowed."
+ write (iout,*) "Specify a smaller r_cut_int and resubmit"
+ call flush(iout)
+ endif
+ write (*,*) "Processor:",me,": Too many SCp interactions",
+ & g_ilist_scp_long+g_ilist_scp_short," only",
+ & 2*maxres*maxint_res," allowed."
+ write (*,*) "Specify a smaller r_cut_int and resubmit"
+ call MPI_Abort(MPI_COMM_WORLD,ierr)
+ endif
+c write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_scp_long,1,MPI_INTEGER,
+ & i_ilist_scp_long,1,MPI_INTEGER,king,FG_COMM,IERR)
+ call MPI_Gather(ilist_scp_short,1,MPI_INTEGER,
+ & i_ilist_scp_short,1,MPI_INTEGER,king,FG_COMM,IERR)
+c write (iout,*) "SCp after gather ierr",ierr
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_scp_long(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistscpi_long,ilist_scp_long,MPI_INTEGER,
+ & newcontlistscpi_long,i_ilist_scp_long,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp after gatherv ierr",ierr
+ call MPI_Gatherv(contlistscpj_long,ilist_scp_long,MPI_INTEGER,
+ & newcontlistscpj_long,i_ilist_scp_long,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp after gatherv ierr",ierr
+ call MPI_Bcast(g_ilist_scp_long,1,MPI_INT,king,FG_COMM,IERR)
+c write (iout,*) "SCp after bcast ierr",ierr
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistscpi_long,g_ilist_scp_long,MPI_INT,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp after bcast ierr",ierr
+ call MPI_Bcast(newcontlistscpj_long,g_ilist_scp_long,MPI_INT,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp bcast reduce ierr",ierr
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_scp_short(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistscpi_short,ilist_scp_short,MPI_INTEGER,
+ & newcontlistscpi_short,i_ilist_scp_short,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp after gatherv ierr",ierr
+ call MPI_Gatherv(contlistscpj_short,ilist_scp_short,MPI_INTEGER,
+ & newcontlistscpj_short,i_ilist_scp_short,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp after gatherv ierr",ierr
+ call MPI_Bcast(g_ilist_scp_short,1,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistscpi_short,g_ilist_scp_short,MPI_INT,
+ & king,FG_COMM,IERR)
+c write (iout,*) "SCp after bcast ierr",ierr
+ call MPI_Bcast(newcontlistscpj_short,g_ilist_scp_short,MPI_INT,
+ & king,FG_COMM,IERR)
+ else
+#endif
+ g_ilist_scp_long=ilist_scp_long
+
+ do i=1,ilist_scp_long
+ newcontlistscpi_long(i)=contlistscpi_long(i)
+ newcontlistscpj_long(i)=contlistscpj_long(i)
+ enddo
+ g_ilist_scp_short=ilist_scp_short
+
+ do i=1,ilist_scp_short
+ newcontlistscpi_short(i)=contlistscpi_short(i)
+ newcontlistscpj_short(i)=contlistscpj_short(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+ if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+ &then
+ write (iout,'(a30,i10,a,i4)')
+ & "Number of long-range SC-p interactions",
+ & g_ilist_scp_long," per residue on average",g_ilist_scp_long/nres
+ write (iout,'(a30,i10,a,i4)')
+ & "Number of short-range SC-p interactions",
+ &g_ilist_scp_short," per residue on average",g_ilist_scp_short/nres
+ endif
+#ifdef DEBUG
+ write (iout,*) "make_SCp_inter_list: after GATHERV long-range",
+ & g_ilist_scp_long
+ do i=1,g_ilist_scp_long
+ write (iout,*) i,newcontlistscpi_long(i),newcontlistscpj_long(i)
+ enddo
+ write (iout,*) "make_SCp_inter_list: after GATHERV short-range",
+ & g_ilist_scp_short
+ do i=1,g_ilist_scp_short
+ write (iout,*) i,newcontlistscpi_short(i),
+ & newcontlistscpj_short(i)
+ enddo
+#endif
+ call int_bounds(g_ilist_scp_long,g_listscp_start_long,
+ & g_listscp_end_long)
+ call int_bounds(g_ilist_scp_short,g_listscp_start_short,
+ & g_listscp_end_short)
+ if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+ &then
+ write (iout,*) "g_listscp_start",g_listscp_start_long,
+ & "g_listscp_end",g_listscp_end_long
+ write (iout,*)"g_listscp_start_short",g_listscp_start_short,
+ & "g_listscp_end_short",g_listscp_end_short
+ endif
return
end
!-----------------------------------------------------------------------------
- subroutine make_pp_vdw_inter_list
+ subroutine make_pp_vdw_inter_list_RESPA
implicit none
include "DIMENSIONS"
#ifdef MPI
double precision xmedj,ymedj,zmedj
double precision dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,
& xmedi,ymedi,zmedi
- double precision dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,
- & dx_normj,dy_normj,dz_normj
- integer contlistpp_vdwi(maxint_res*maxres),
- & contlistpp_vdwj(maxint_res*maxres)
-! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+ double precision dxj,dyj,dzj
+ integer contlistpp_vdwi_short(maxint_res*maxres),
+ & contlistpp_vdwj_short(maxint_res*maxres)
integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
- & ilist_pp_vdw,g_ilist_pp_vdw
- integer displ(0:max_fg_procs),i_ilist_pp_vdw(0:max_fg_procs),ierr
+ & ilist_pp_vdw_short,g_ilist_pp_vdw_short
+ integer displ(0:max_fg_procs),
+ & i_ilist_pp_vdw_short(0:max_fg_procs),ierr
! print *,"START make_SC"
+ double precision boxshift
+ double precision d_scale,r_respa_buf
#ifdef DEBUG
write (iout,*) "make_pp_vdw_inter_list"
#endif
- ilist_pp_vdw=0
+ ilist_pp_vdw_short=0
r_buff_list=5.0
+ r_respa_buf=rlamb
do i=iatel_s_vdw,iatel_e_vdw
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
- dx_normi=dc_norm(1,i)
- dy_normi=dc_norm(2,i)
- dz_normi=dc_norm(3,i)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=dmod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=dmod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=dmod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
do j=ielstart_vdw(i),ielend_vdw(i)
! write (iout,*) i,j,itype(i),itype(j)
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
-
! 1,j)
dxj=dc(1,j)
dyj=dc(2,j)
dzj=dc(3,j)
- dx_normj=dc_norm(1,j)
- dy_normj=dc_norm(2,j)
- dz_normj=dc_norm(3,j)
! xj=c(1,j)+0.5D0*dxj-xmedi
! yj=c(2,j)+0.5D0*dyj-ymedi
! zj=c(3,j)+0.5D0*dzj-zmedi
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
-
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- endif
- enddo
- enddo
- enddo
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
+ dist_init=dsqrt(xj*xj+yj*yj+zj*zj)
- if (dsqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+ if (dist_init.le.(r_cut_int+r_buff_list)) then
+ d_scale=dist_init/rpp(itel(i),itel(j))
+ if (d_scale.le.r_cut_respa+r_respa_buf) then
! Here the list is created
- ilist_pp_vdw=ilist_pp_vdw+1
+ ilist_pp_vdw_short=ilist_pp_vdw_short+1
! this can be substituted by cantor and anti-cantor
- contlistpp_vdwi(ilist_pp_vdw)=i
- contlistpp_vdwj(ilist_pp_vdw)=j
+ contlistpp_vdwi_short(ilist_pp_vdw_short)=i
+ contlistpp_vdwj_short(ilist_pp_vdw_short)=j
+ endif
endif
- enddo
- enddo
+ enddo
+ enddo
! enddo
#ifdef MPI
#ifdef DEBUG
- write (iout,*) "before MPIREDUCE",ilist_pp_vdw
- do i=1,ilist_pp_vdw
- write (iout,*) i,contlistpp_vdwi(i),contlistpp_vdwj(i)
+ write (iout,*) "before MPIREDUCE longrange",ilist_pp_vdw_long
+ do i=1,ilist_pp_vdw_long
+ write (iout,*) i,contlistpp_vdwi_long(i),contlistpp_vdwj_long(i)
+ enddo
+ write (iout,*) "before MPIREDUCE shortrange",ilist_pp_vdw_short
+ do i=1,ilist_pp_vdw_short
+ write (iout,*) i,contlistpp_vdwi_short(i),
+ & contlistpp_vdwj_short(i)
enddo
#endif
if (nfgtasks.gt.1)then
- call MPI_Reduce(ilist_pp_vdw,g_ilist_pp_vdw,1,
+ call MPI_Reduce(ilist_pp_vdw_short,g_ilist_pp_vdw_short,1,
& MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
- if (fg_rank.eq.0.and.g_ilist_pp_vdw.gt.maxres*maxint_res) then
+ if (fg_rank.eq.0.and.g_ilist_pp_vdw_short.gt.maxres*maxint_res)
+ & then
if ((me.eq.king.or.out1file).and.energy_dec) then
write (iout,*) "Too many pp VDW interactions",
- & g_ilist_pp_vdw," only",maxres*maxint_res," allowed."
+ & g_ilist_pp_vdw_short," only",maxres*maxint_res," allowed."
write (iout,*) "Specify a smaller r_cut_int and resubmit"
call flush(iout)
endif
write (*,*) "Processor:",me,": Too many pp VDW interactions",
- & g_ilist_pp_vdw," only",maxres*maxint_res," allowed."
+ & g_ilist_pp_vdw_short," only",maxres*maxint_res," allowed."
write (8,*) "Specify a smaller r_cut_int and resubmit"
call MPI_Abort(MPI_COMM_WORLD,ierr)
endif
! write(iout,*) "before bcast",g_ilist_sc
- call MPI_Gather(ilist_pp_vdw,1,MPI_INTEGER,
- & i_ilist_pp_vdw,1,MPI_INTEGER,king,FG_COMM,IERR)
+ call MPI_Gather(ilist_pp_vdw_short,1,MPI_INTEGER,
+ & i_ilist_pp_vdw_short,1,MPI_INTEGER,king,FG_COMM,IERR)
displ(0)=0
do i=1,nfgtasks-1,1
- displ(i)=i_ilist_pp_vdw(i-1)+displ(i-1)
+ displ(i)=i_ilist_pp_vdw_short(i-1)+displ(i-1)
enddo
! write(iout,*) "before gather",displ(0),displ(1)
- call MPI_Gatherv(contlistpp_vdwi,ilist_pp_vdw,MPI_INTEGER,
- & newcontlistpp_vdwi,i_ilist_pp_vdw,displ,MPI_INTEGER,
- & king,FG_COMM,IERR)
- call MPI_Gatherv(contlistpp_vdwj,ilist_pp_vdw,MPI_INTEGER,
- & newcontlistpp_vdwj,i_ilist_pp_vdw,displ,MPI_INTEGER,
- & king,FG_COMM,IERR)
- call MPI_Bcast(g_ilist_pp_vdw,1,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistpp_vdwi_short,ilist_pp_vdw_short,
+ & MPI_INTEGER,newcontlistpp_vdwi_short,i_ilist_pp_vdw_short,displ,
+ & MPI_INTEGER,king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistpp_vdwj_short,ilist_pp_vdw_short,
+ & MPI_INTEGER,newcontlistpp_vdwj_short,i_ilist_pp_vdw_short,displ,
+ & MPI_INTEGER,king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_pp_vdw_short,1,MPI_INT,king,FG_COMM,IERR)
! write(iout,*) "before bcast",g_ilist_sc
! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
- call MPI_Bcast(newcontlistpp_vdwi,g_ilist_pp_vdw,MPI_INT,king,
- & FG_COMM,IERR)
- call MPI_Bcast(newcontlistpp_vdwj,g_ilist_pp_vdw,MPI_INT,king,
- & FG_COMM,IERR)
-
+ call MPI_Bcast(newcontlistpp_vdwi_short,g_ilist_pp_vdw_short,
+ & MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistpp_vdwj_short,g_ilist_pp_vdw_short,
+ & MPI_INT,king,FG_COMM,IERR)
! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
-
- else
+ else
#endif
- g_ilist_pp_vdw=ilist_pp_vdw
+ g_ilist_pp_vdw_short=ilist_pp_vdw_short
- do i=1,ilist_pp_vdw
- newcontlistpp_vdwi(i)=contlistpp_vdwi(i)
- newcontlistpp_vdwj(i)=contlistpp_vdwj(i)
+ do i=1,ilist_pp_vdw_short
+ newcontlistpp_vdwi_short(i)=contlistpp_vdwi_short(i)
+ newcontlistpp_vdwj_short(i)=contlistpp_vdwj_short(i)
enddo
#ifdef MPI
- endif
+ endif
#endif
- call int_bounds(g_ilist_pp_vdw,g_listpp_vdw_start,
- & g_listpp_vdw_end)
- if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
- &write (iout,'(a30,i10,a,i4)') "Number of p-p VDW interactions",
- & g_ilist_pp_vdw," per residue on average",g_ilist_pp_vdw/nres
+ if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+ &then
+ write (iout,*) "Number of short-range p-p VDW interactions",
+ & g_ilist_pp_vdw_short," per residue on average",
+ & g_ilist_pp_vdw_short/nres
+ endif
#ifdef DEBUG
- write (iout,*) "g_listpp_vdw_start",g_listpp_vdw_start,
- & "g_listpp_vdw_end",g_listpp_vdw_end
+ write (iout,*) "Short-range pp_vdw"
write (iout,*) "make_pp_vdw_inter_list: after GATHERV",
- & g_ilist_pp_vdw
- do i=1,g_ilist_pp_vdw
- write (iout,*) i,newcontlistpp_vdwi(i),newcontlistpp_vdwj(i)
+ & g_ilist_pp_vdw_short
+ do i=1,g_ilist_pp_vdw_short
+ write (iout,*) i,newcontlistpp_vdwi_short(i),
+ & newcontlistpp_vdwj_short(i)
enddo
#endif
+ call int_bounds(g_ilist_pp_vdw_short,g_listpp_vdw_start_short,
+ & g_listpp_vdw_end_short)
+ if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+ &then
+ write (iout,*)"g_listpp_vdw_start_short",
+ & g_listpp_vdw_start_short,
+ & "g_listpp_vdw_end_short",g_listpp_vdw_end_short
+ endif
return
end
!-----------------------------------------------------------------------------
! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
- else
+ else
#endif
g_ilist_pp=ilist_pp
do i=1,ilist_pp
- newcontlistppi(i)=contlistppi(i)
- newcontlistppj(i)=contlistppj(i)
+ newcontlistppi(i)=contlistppi(i)
+ newcontlistppj(i)=contlistppj(i)
enddo
#ifdef MPI
- endif
+ endif
#endif
- call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
& write (iout,'(a30,i10,a,i4)') "Number of p-p interactions",
& g_ilist_pp," per residue on average",g_ilist_pp/nres
write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
enddo
#endif
+ call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+ if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+ &then
+ write (iout,*) "g_listpp_start",g_listpp_start,
+ & "g_listpp_end",g_listpp_end
+ endif
return
end
#endif
include 'COMMON.TIME1'
double precision z(maxres6),d_a_tmp(maxres6)
- double precision edum(0:n_ene),time_order(0:10)
+ double precision edum(0:n_ene),time_order(0:11)
c double precision Gcopy(maxres2,maxres2)
c common /przechowalnia/ Gcopy
integer icall /0/
- integer i,j,iorder
+ integer i,j,iorder,ioverlap(maxres),ioverlap_last
C Workers wait for variables and NF, and NFL from the boss
iorder=0
do while (iorder.ge.0)
#endif
else if (iorder.eq.10) then
call setup_fricmat
+ else if (iorder.eq.11) then
+ call overlap_sc_list(ioverlap,ioverlap_last,.false.)
endif
enddo
write (*,*) 'Processor',fg_rank,' CG group',kolor,
--- /dev/null
+ subroutine orig_frame_chain(istart)
+C
+C Define the origin and orientation of the coordinate system starting
+C at residue istart and locate sites istart+1 and istart+2. The
+C coordinates of site istart and the respective dc and dc_norm must be
+C pre-defined
+C
+ implicit none
+ integer i,j,istart,ichain
+ double precision cost,sint,cosg,sing,aux
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.LOCAL'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ cost=dc_norm(1,istart)
+ aux=dsqrt(dc_norm(2,istart)**2+dc_norm(3,istart)**2)
+ cosg=dc_norm(2,istart)/aux
+ sing=dc_norm(3,istart)/aux
+ cost=dcos(theta(istart+2))
+ sint=dsin(theta(istart+2))
+ prod(1,1,istart)=cost
+ prod(1,2,istart)=-sint
+ prod(1,3,istart)=0.0d0
+ prod(2,1,istart)=sint*cosg
+ prod(2,2,istart)=cost*cosg
+ prod(2,3,istart)=-sing
+ prod(3,1,istart)=sint*sing
+ prod(3,2,istart)=cost*sing
+ prod(3,3,istart)=cosg
+ aux=prod(1,1,istart)*(prod(2,2,istart)*prod(3,3,istart)
+ & -prod(3,2,istart)*prod(2,3,istart))
+ & -prod(1,2,istart)*(prod(2,1,istart)*prod(3,3,istart)
+ & -prod(3,1,istart)*prod(2,3,istart))
+ & +prod(1,3,istart)*(prod(2,1,istart)*prod(3,2,istart)
+ & -prod(3,1,istart)*prod(2,2,istart))
+c write (iout,*) "orig_frame_chain prod",istart
+c do i=1,3
+c write(iout,'(i5,3f10.5)') i,(prod(i,j,istart),j=1,3)
+c enddo
+c write (iout,*) "orig_frame_chain: prod",istart," determinant",aux
+ t(1,1,istart)=cost
+ t(1,2,istart)=-sint
+ t(1,3,istart)= 0.0D0
+ t(2,1,istart)=sint
+ t(2,2,istart)=cost
+ t(2,3,istart)= 0.0D0
+ t(3,1,istart)= 0.0D0
+ t(3,2,istart)= 0.0D0
+ t(3,3,istart)= 1.0D0
+ r(1,1,istart)= 1.0D0
+ r(1,2,istart)= 0.0D0
+ r(1,3,istart)= 0.0D0
+ r(2,1,istart)= 0.0D0
+ r(2,2,istart)= 1.0D0
+ r(2,3,istart)= 0.0D0
+ r(3,1,istart)= 0.0D0
+ r(3,2,istart)= 0.0D0
+ r(3,3,istart)= 1.0D0
+ do i=1,3
+ do j=1,3
+ rt(i,j,istart)=t(i,j,istart)
+ enddo
+ enddo
+ call matmult(prod(1,1,istart),rt(1,1,istart),prod(1,1,istart+1))
+c aux=prod(1,1,istart+1)*(prod(2,2,istart+1)*prod(3,3,istart+1)
+c & -prod(3,2,istart+1)*prod(2,3,istart+1))
+c & -prod(1,2,istart+1)*(prod(2,1,istart+1)*prod(3,3,istart+1)
+c & -prod(3,1,istart+1)*prod(2,3,istart+1))
+c & +prod(1,3,istart+1)*(prod(2,1,istart+1)*prod(3,2,istart+1)
+c & -prod(3,1,istart+1)*prod(2,2,istart+1))
+c write (iout,*) "orig_frame_chain prod",istart+1
+c do i=1,3
+c write(iout,'(i5,3f10.5)') i,(prod(i,j,istart+1),j=1,3)
+c enddo
+c write (iout,*)"orig_frame_chain: prod",istart+1," determinant",aux
+ do j=1,3
+ dc_norm(j,istart+1)=prod(j,1,istart+1)
+ dc(j,istart+1)=vbld(istart+2)*prod(j,1,istart+1)
+ c(j,istart+2)=c(j,istart+1)+dc(j,istart+1)
+ enddo
+ call locate_side_chain(istart+1)
+ return
+ end
iterter(ires_old-1)=1
itype(ires_old)=ntyp1
iterter(ires_old)=1
- ishift1=ishift1+1
+c ishift1=ishift1+1
ibeg=2
- write (iout,*) "Chain ended",ires,ishift,ires_old
+ write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
read (card(18:20),'(a3)') res
c write (iout,*) "ires",ires,ires-ishift+ishift1,
c & " ires_old",ires_old
-c write (iout,*) "ishift",ishift," ishift1",ishift1
-c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
if (ires-ishift+ishift1.ne.ires_old) then
! Calculate the CM of the preceding residue.
! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
sccalc=.true.
endif
! Start new residue.
+c write (iout,*) "ibeg",ibeg
if (res.eq.'Cl-' .or. res.eq.'Na+') then
ires=ires_old
cycle
else if (ibeg.eq.2) then
! Start a new chain
ishift=-ires_old+ires-1 !!!!!
- ishift1=ishift1-1 !!!!!
-c write (iout,*) "New chain started",ires,ishift,ishift1,"!"
+c ishift1=ishift1-1 !!!!!
+c write (iout,*) "New chain started",ires,ires_old,ishift,
+c & ishift1
ires=ires-ishift+ishift1
+ write (iout,*) "New chain started ires",ires
ires_old=ires
+c ires=ires_old+1
ibeg=0
else
ishift=ishift-(ires-ishift+ishift1-ires_old-1)
if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
& res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-! write (iout,*) "backbone ",atom
+ read(card(61:66),*) bfac(ires)
+c write (iout,*) "backbone ",atom
+c write (iout,*) ires,res,(c(j,ires),j=1,3)
#ifdef DEBUG
write (iout,'(i6,i3,2x,a,3f8.3)')
& ires,itype(ires),res,(c(j,ires),j=1,3)
endif
enddo
10 if(me.eq.king.or..not.out1file)
- & write (iout,'(a,i5)') ' Nres: ',ires
+ & write (iout,'(a,i7)') ' Nres: ',ires
c write (iout,*) "iii",iii
C Calculate dummy residue coordinates inside the "chain" of a multichain
C system
e2(3)=0.0d0
endif
do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
enddo
else !unres_pdb
do j=1,3
write (iout,'(/a)')
& "Cartesian coordinates of the reference structure"
write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
- & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
do ires=1,nres
- write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')
+ write (iout,'(a3,1x,i6,3f8.3,5x,3f8.3)')
& restyp(itype(ires)),ires,(c(j,ires),j=1,3),
& (c(j,ires+nres),j=1,3)
enddo
character*3 seq,res
character*5 atom
character*80 card
- double precision sccor(3,20)
+ double precision sccor(3,50)
integer rescode,iterter(maxres)
+ logical zero
do i=1,maxres
iterter(i)=0
enddo
endif
enddo
10 if(me.eq.king.or..not.out1file)
- & write (iout,'(a,i5)') ' Nres: ',ires
+ & write (iout,'(a,i7)') ' Nres: ',ires
C Calculate dummy residue coordinates inside the "chain" of a multichain
C system
nres=ires
e2(3)=0.0d0
endif
do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
enddo
else !unres_pdb
do j=1,3
e2(3)=0.0d0
endif
do j=1,3
- c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+ c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
enddo
else
do j=1,3
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-1.9d0*e2(j)
+ c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
enddo
else
do j=1,3
& (c(j,ires+nres),j=1,3)
enddo
endif
+ zero=.false.
+ do ires=1,nres
+ zero=zero.or.itype(ires).eq.0
+ enddo
+ if (zero) then
+ write (iout,'(2a)') "Gaps in PDB coordinates detected;",
+ & " look for ZERO in the control output above."
+ write (iout,'(2a)') "Repair the PDB file using MODELLER",
+ & " or other softwared and resubmit."
+ call flush(iout)
+ stop
+ endif
C Calculate internal coordinates.
call int_from_cart(.true.,out_template_coord)
call sc_loc_geom(.false.)
thetaref(i)=theta(i)
phiref(i)=phi(i)
enddo
+ dc(:,0)=c(:,1)
do i=1,nres-1
do j=1,3
dc(j,i)=c(j,i+1)-c(j,i)
out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
call readi(controlcard,'SYM',symetr,1)
+ call readi(controlcard,'PERMUT',npermut,1)
call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
c endif
call readi(controlcard,'NZ_START',nz_start,0)
call readi(controlcard,'NZ_END',nz_end,0)
+ call readi(controlcard,'NRAN_START',nran_start,0)
+ write (iout,*) "nran_start",nran_start
c call readi(controlcard,'IZ_SC',iz_sc,0)
timlim=60.0D0*timlim
safety = 60.0d0*safety
call readi(controlcard,"NSTEP",n_timestep,1000000)
call readi(controlcard,"NTWE",ntwe,100)
call readi(controlcard,"NTWX",ntwx,1000)
+ call readi(controlcard,"REST_FREQ",irest_freq,1000)
call reada(controlcard,"DT",d_time,1.0d-1)
call reada(controlcard,"DVMAX",dvmax,2.0d1)
call reada(controlcard,"DAMAX",damax,1.0d1)
& " A"
write(iout,'(a60,i5)')"Frequency of updating interaction list",
& imatupdate
+ write(iout,'(a60,i5)')"Restart writing frequency",irest_freq
if (RESPA) then
write (iout,'(2a,i4,a)')
& "A-MTS algorithm used; initial time step for fast-varying",
integer ilen
external ilen
integer iperm,tperm
- integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp
+ integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
double precision sumv
C
C Read PDB structure if applicable
& chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
enddo
call chain_symmetry(nchain,nres,itype,chain_border,
- & chain_length,npermchain,tabpermchain)
+ & chain_length,npermchain,tabpermchain,nchain_group,nequiv,
+ & iequiv,mapchain)
c do i=1,nres
c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
c & ii=1,npermchain)
c enddo
+#ifdef DEBUG
write(iout,*) "residue permutations"
do i=1,nres
write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
enddo
call flush(iout)
+#endif
if (itype(1).eq.ntyp1) nnt=2
if (itype(nres).eq.ntyp1) nct=nct-1
write (iout,*) "nnt",nnt," nct",nct
call reada(weightcard,"WDIHC",wdihc,0.591D0)
write (iout,*) "Weight of dihedral angle restraints",wdihc
read(inp,'(9x,3f7.3)')
+c & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
& (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
write (iout,*) "The secprob array"
do i=nnt,nct
write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
enddo
ndih_constr=0
+ iconstr_dih=0
do i=nnt+3,nct
if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
& .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
ndih_constr=ndih_constr+1
idih_constr(ndih_constr)=i
+ iconstr_dih(i)=ndih_constr
sumv=0.0d0
do j=1,3
vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
sumv=sumv+vpsipred(j,ndih_constr)
enddo
- do j=1,3
- vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
- enddo
+ if (sumv.gt.0.0d0) then
+ do j=1,3
+ vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+ enddo
+ else
+ vpsipred(1,ndih_constr)=1.0d0
+ vpsipred(2,ndih_constr)=0.0d0
+ vpsipred(3,ndih_constr)=0.0d0
+ endif
phibound(1,ndih_constr)=phihel*deg2rad
phibound(2,ndih_constr)=phibet*deg2rad
sdihed(1,ndih_constr)=sigmahel*deg2rad
endif
write (iout,*) "calling read_saxs_consrtr",nsaxs
if (nsaxs.gt.0) call read_saxs_constr
-
+c write (iout,*) "After read_saxs_constr"
+c call flush(iout)
if (constr_homology.gt.0) then
+c write (iout,*) "Calling read_constr_homology"
+c call flush(iout)
call read_constr_homology
if (indpdb.gt.0 .or. pdbref) then
do i=1,2*nres
335 continue
unres_pdb=.false.
nres_temp=nres
- call readpdb
+c call readpdb
+ call readpdb_template(nmodel_start+1)
close(ipdbin)
if (nres.ge.nres_temp) then
nmodel_start=nmodel_start+1
enddo
enddo
else
-c itemp=nres
-c nres=nres_temp
-c call gen_rand_conf(itemp,*115)
-c nmodel_start=nmodel_start+1
-c do i=1,2*nres
-c do j=1,3
-c chomo(j,i,nmodel_start)=c(j,i)
-c enddo
-c enddo
-c goto 116
- 115 if (me.eq.king .or. .not. out1file)
- & write (iout,'(a,2i5,1x,a)')
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,'(a,2i7,1x,a)')
& "Different number of residues",nres_temp,nres,
& " model skipped."
endif
- 116 continue
nres=nres_temp
enddo
332 continue
#else
if (me.eq.king .or. .not.out1file)
& open(iout,file=outname,status='unknown')
+#define DEBUG
#ifdef DEBUG
if (fg_rank.gt.0) then
write (liczba,'(i3.3)') myrank/nfgtasks
& status='unknown')
endif
#endif
+#undef DEBUG
if(me.eq.king) then
open(igeom,file=intname,status='unknown',access='append')
open(ipdb,file=pdbname,status='unknown')
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#else
- write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
enddo
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#else
- write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
enddo
ipeak(2,npeak)=i
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
& nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
& ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
& dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
& fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
#else
- write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
& nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
& ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
& dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
irestr_type(nhpb)=11
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
& dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
#else
- write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
& dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
#endif
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ & write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
& dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
& irestr_type(nhpb)
#else
- write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
& dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
& irestr_type(nhpb)
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ & write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#else
- write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#endif
endif
dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
#ifdef MPI
if (.not.out1file .or. me.eq.king) then
- write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
& dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
& irestr_type(nhpb)
endif
#else
- write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
& dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
& irestr_type(nhpb)
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
- & ik,iistart,nres_temp
+ integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp,
+ & irec,ik,iistart,nres_temp
integer ilen
external ilen
logical liiflag,lfirst
call card_concat(controlcard)
read(controlcard,*) (waga_homology(i),i=1,homol_nset)
if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
- write(iout,*) "iset homology_weight "
+c write(iout,*) "iset homology_weight "
do i=1,homol_nset
write(iout,*) i,waga_homology(i)
enddo
cd write (iout,*) "nnt",nnt," nct",nct
cd call flush(iout)
+ if (read_homol_frag) then
+ call read_klapaucjusz
+ else
lim_odl=0
lim_dih=0
c
c write(iout,*) 'nnt=',nnt,'nct=',nct
c
- do i = nnt,nct
- do k=1,constr_homology
- idomain(k,i)=0
- enddo
- enddo
-
- ii=0
- do i = nnt,nct-2
- do j=i+2,nct
- ii=ii+1
- ii_in_use(ii)=0
- enddo
- enddo
-
- if (read_homol_frag) then
- call read_klapaucjusz
- else
+c do i = nnt,nct
+c do k=1,constr_homology
+c idomain(k,i)=0
+c enddo
+c enddo
+ idomain=0
+c ii=0
+c do i = nnt,nct-2
+c do j=i+2,nct
+c ii=ii+1
+c ii_in_use(ii)=0
+c enddo
+c enddo
+ ii_in_use=0
do k=1,constr_homology
read(inp,'(a)') pdbfile
endif
sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
else
- ii=ii+1
- l_homo(k,ii)=.false.
+c ii=ii+1
+c l_homo(k,ii)=.false.
endif
enddo
enddo
if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write (iout,*) "Distance restraints from templates"
do ii=1,lim_odl
- write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
+ write(iout,'(3i7,100(2f8.2,1x,l1,4x))')
& ii,ires_homo(ii),jres_homo(ii),
& (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
& ki=1,constr_homology)
enddo
write (iout,*) "Dihedral angle restraints from templates"
do i=nnt+3,nct
- write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+ write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
& (rad2deg*dih(ki,i),
& rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
enddo
write (iout,*) "Virtual-bond angle restraints from templates"
do i=nnt+2,nct
- write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+ write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
& (rad2deg*thetatpl(ki,i),
& rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
enddo
write (iout,*) "SC restraints from templates"
do i=nnt,nct
- write(iout,'(i5,100(4f8.2,4x))') i,
+ write(iout,'(i7,100(4f8.2,4x))') i,
& (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
& 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
enddo
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
- & ik,ll,ii,kk,iistart,iishift,lim_xx
+ integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
+ & idomain_tmp,
+ & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
+ & i01,i10,nnt_chain,nct_chain
+ integer itype_temp(maxres)
double precision distal
logical lprn /.true./
integer nres_temp
integer ilen
external ilen
- logical liiflag
+ logical liiflag,lfirst
c
c
double precision rescore_tmp,x12,y12,z12,rescore2_tmp
c
include 'COMMON.VAR'
c
+c write (iout,*) "READ_KLAPAUCJUSZ"
+c print *,"READ_KLAPAUCJUSZ"
+c call flush(iout)
call getenv("FRAGFILE",fragfile)
+ write (iout,*) "Opening", fragfile
+ call flush(iout)
open(ientin,file=fragfile,status="old",err=10)
- read(ientin,*) constr_homology,nclust
- nmodel_start=constr_homology
- l_homo = .false.
+c write (iout,*) " opened"
+c call flush(iout)
+
sigma_theta=0.0
sigma_d=0.0
sigma_dih=0.0
+ l_homo = .false.
+
+ nres_temp=nres
+ itype_temp=itype
+ ii=0
+ lim_odl=0
+
+c write (iout,*) "Entering loop"
+c call flush(iout)
+
+ DO IGR = 1,NCHAIN_GROUP
+
+c write (iout,*) "igr",igr
+ call flush(iout)
+ read(ientin,*) constr_homology,nclust
+
+ if (start_from_model) then
+ nmodel_start=constr_homology
+ else
+ nmodel_start=0
+ endif
+
+ ii_old=lim_odl
+
+ ichain=iequiv(1,igr)
+ nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
+ nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
+c write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
+
c Read pdb files
do k=1,constr_homology
read(ientin,'(a)') pdbfile
stop
34 continue
unres_pdb=.false.
- nres_temp=nres
+c nres_temp=nres
call readpdb_template(k)
nres_chomo(k)=nres
- nres=nres_temp
+c nres=nres_temp
do i=1,nres
rescore(k,i)=0.2d0
rescore2(k,i)=1.0d0
do ll = 1,ninclust(l)
k = inclust(ll,l)
+c write (iout,*) "l",l," ll",ll," k",k
do i=1,nres
idomain(k,i)=0
enddo
c
c ... --> odl(k,ii)
C Copy the coordinates from reference coordinates (?)
- nres_temp=nres
+c nres_temp=nres
nres=nres_chomo(k)
do i=1,2*nres
do j=1,3
thetaref(i)=theta(i)
phiref(i)=phi(i)
enddo
- nres=nres_temp
+c nres=nres_temp
if (waga_dist.ne.0.0d0) then
- ii=0
- do i = nnt,nct-2
- do j=i+2,nct
+ ii=ii_old
+c do i = nnt,nct-2
+ do i = nnt_chain,nct_chain-2
+c do j=i+2,nct
+ do j=i+2,nct_chain
x12=c(1,i)-c(1,j)
y12=c(2,i)-c(2,j)
c write (iout,*) "k",k
c write (iout,*) "i",i," j",j," constr_homology",
-c & constr_homology
- ires_homo(ii)=i
- jres_homo(ii)=j
+c & constr_homology
+ ires_homo(ii)=i+chain_border1(1,igr)-1
+ jres_homo(ii)=j+chain_border1(1,igr)-1
odl(k,ii)=distal
if (read2sigma) then
sigma_odl(k,ii)=0
c Theta, dihedral and SC retraints
c
if (waga_angle.gt.0.0d0) then
- do i = nnt+3,nct
+ do i = nnt_chain+3,nct_chain
+ iii=i+chain_border1(1,igr)-1
if (idomain(k,i).eq.0) then
c sigma_dih(k,i)=0.0
cycle
endif
- dih(k,i)=phiref(i)
- sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ dih(k,iii)=phiref(i)
+ sigma_dih(k,iii)=
+ & (rescore(k,i)+rescore(k,i-1)+
& rescore(k,i-2)+rescore(k,i-3))/4.0
c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
c & " sigma_dihed",sigma_dih(k,i)
- if (sigma_dih(k,i).ne.0)
- & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+ if (sigma_dih(k,iii).ne.0)
+ & sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
enddo
- lim_dih=nct-nnt-2
+c lim_dih=nct-nnt-2
endif
if (waga_theta.gt.0.0d0) then
- do i = nnt+2,nct
+ do i = nnt_chain+2,nct_chain
+ iii=i+chain_border1(1,igr)-1
if (idomain(k,i).eq.0) then
c sigma_theta(k,i)=0.0
cycle
endif
- thetatpl(k,i)=thetaref(i)
- sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ thetatpl(k,iii)=thetaref(i)
+ sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+
& rescore(k,i-2))/3.0
- if (sigma_theta(k,i).ne.0)
- & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+ if (sigma_theta(k,iii).ne.0)
+ & sigma_theta(k,iii)=1.0d0/
+ & (sigma_theta(k,iii)*sigma_theta(k,iii))
enddo
endif
if (waga_d.gt.0.0d0) then
- do i = nnt,nct
+ do i = nnt_chain,nct_chain
+ iii=i+chain_border1(1,igr)-1
if (itype(i).eq.10) cycle
if (idomain(k,i).eq.0 ) then
c sigma_d(k,i)=0.0
cycle
endif
- xxtpl(k,i)=xxref(i)
- yytpl(k,i)=yyref(i)
- zztpl(k,i)=zzref(i)
- sigma_d(k,i)=rescore(k,i)
- if (sigma_d(k,i).ne.0)
- & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
- if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
+ xxtpl(k,iii)=xxref(i)
+ yytpl(k,iii)=yyref(i)
+ zztpl(k,iii)=zzref(i)
+ sigma_d(k,iii)=rescore(k,i)
+ if (sigma_d(k,iii).ne.0)
+ & sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
+c if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
enddo
endif
enddo ! l
c remove distance restraints not used in any model from the list
c shift data in all arrays
c
+c write (iout,*) "ii_old",ii_old
if (waga_dist.ne.0.0d0) then
- ii=0
+#ifdef DEBUG
+ write (iout,*) "Distance restraints from templates"
+ do iii=1,lim_odl
+ write(iout,'(4i5,100(2f8.2,1x,l1,4x))')
+ & iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
+ & (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
+ & ki=1,constr_homology)
+ enddo
+#endif
+ ii=ii_old
liiflag=.true.
- do i=nnt,nct-2
- do j=i+2,nct
+ lfirst=.true.
+ do i=nnt_chain,nct_chain-2
+ do j=i+2,nct_chain
ii=ii+1
- if (ii_in_use(ii).eq.0.and.liiflag) then
+c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c & .and. distal.le.dist2_cut ) then
+c write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
+c call flush(iout)
+ if (ii_in_use(ii).eq.0.and.liiflag.or.
+ & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
liiflag=.false.
- iistart=ii
+ i10=ii
+ if (lfirst) then
+ lfirst=.false.
+ iistart=ii
+ else
+ if(i10.eq.lim_odl) i10=i10+1
+ do ki=0,i10-i01-1
+ ires_homo(iistart+ki)=ires_homo(ki+i01)
+ jres_homo(iistart+ki)=jres_homo(ki+i01)
+ ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+ do k=1,constr_homology
+ odl(k,iistart+ki)=odl(k,ki+i01)
+ sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+ l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+ enddo
+ enddo
+ iistart=iistart+i10-i01
+ endif
endif
- if (ii_in_use(ii).ne.0.and..not.liiflag.or.
- & .not.liiflag.and.ii.eq.lim_odl) then
- if (ii.eq.lim_odl) then
- iishift=ii-iistart+1
- else
- iishift=ii-iistart
- endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag) then
+ i01=ii
liiflag=.true.
- do ki=iistart,lim_odl-iishift
- ires_homo(ki)=ires_homo(ki+iishift)
- jres_homo(ki)=jres_homo(ki+iishift)
- ii_in_use(ki)=ii_in_use(ki+iishift)
- do k=1,constr_homology
- odl(k,ki)=odl(k,ki+iishift)
- sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
- l_homo(k,ki)=l_homo(k,ki+iishift)
- enddo
- enddo
- ii=ii-iishift
- lim_odl=lim_odl-iishift
endif
enddo
enddo
+ lim_odl=iistart-1
endif
+ lll=lim_odl-ii_old
+
+ do i=2,nequiv(igr)
+
+ ichain=iequiv(i,igr)
+
+ do j=nnt_chain,nct_chain
+ jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
+ do k=1,constr_homology
+ dih(k,jj)=dih(k,j)
+ sigma_dih(k,jj)=sigma_dih(k,j)
+ thetatpl(k,jj)=thetatpl(k,j)
+ sigma_theta(k,jj)=sigma_theta(k,j)
+ xxtpl(k,jj)=xxtpl(k,j)
+ yytpl(k,jj)=yytpl(k,j)
+ zztpl(k,jj)=zztpl(k,j)
+ sigma_d(k,jj)=sigma_d(k,j)
+ enddo
+ enddo
+
+ jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
+c write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
+ do j=ii_old+1,lim_odl
+ ires_homo(j+lll)=ires_homo(j)+jj
+ jres_homo(j+lll)=jres_homo(j)+jj
+ do k=1,constr_homology
+ odl(k,j+lll)=odl(k,j)
+ sigma_odl(k,j+lll)=sigma_odl(k,j)
+ l_homo(k,j+lll)=l_homo(k,j)
+ enddo
+ enddo
+
+ ii_old=ii_old+lll
+ lim_odl=lim_odl+lll
+
+ enddo
+
+ ENDDO ! IGR
+
+ if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
+ nres=nres_temp
+ itype=itype_temp
+
return
10 stop "Error in fragment file"
end
c print *,fail
s1=0.0
s2=0.0
- print *,s1,s2
+c print *,s1,s2
do 1 i=1,3
- print *, i2,i3,i4
+c print *, i2,i3,i4
zi=c(i,i2)-c(i,i3)
ui=c(i,i4)-c(i,i3)
- print *,zi,ui
+c print *,zi,ui
s1=s1+zi*zi
s2=s2+ui*ui
z(i)=zi
write(iout,1000) i3,i4,i1
fail=.true.
return
- print *,'two if pass'
+c print *,'two if pass'
4 s1=1.0/s1
s2=1.0/s2
v1=z(2)*u(3)-z(3)*u(2)
subroutine friction_force
implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
include 'DIMENSIONS'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
logical lprn /.false./, checkmode /.false./
#ifdef FIVEDIAG
+#ifdef TIMING
+ include 'COMMON.TIME1'
+ double precision time01
+#endif
c Here accelerations due to friction forces are computed right after forces.
- d_t_work=0.0d0
+ d_t_work(:6*nres)=0.0d0
do j=1,3
v_work(j,1)=d_t(j,0)
v_work(j,nnt)=d_t(j,0)
write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
#endif
c write (iout,*) "chain",i," ind",ind," n",n
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
& DU2fric(iposc),vvec(iposc),rs)
+#ifdef TIMING
+ time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
+#endif
#ifdef DEBUG
write (iout,*) "rs"
write (iout,'(f10.5)') (rs(i),i=1,n)
write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
#endif
#else
- do i=0,MAXRES2
+ do i=0,2*nres
do j=1,3
friction(j,i)=0.0d0
enddo
integer ichain,innt,inct,iposc
#endif
- do i=0,MAXRES2
+ do i=0,2*nres
do j=1,3
stochforc(j,i)=0.0d0
enddo
enddo
endif
#ifdef FIVEDIAG
- DMfric=0.0d0
- DU1fric=0.0d0
- DU2fric=0.0d0
+ DMfric(:2*nres)=0.0d0
+ DU1fric(:2*nres)=0.0d0
+ DU2fric(:2*nres)=0.0d0
ind=1
do ichain=1,nchain
innt=chain_border(1,ichain)
time_fricmatmult=0.0d0
time_fsample=0.0d0
time_SAXS=0.0d0
+ time_list=0.0d0
+ time_evdw=0.0d0
+ time_evdw_short=0.0d0
+ time_evdw_long=0.0d0
+ time_eelec=0.0d0
+ time_eelec_short=0.0d0
+ time_eelec_long=0.0d0
+ time_escp=0.0d0
+ time_escp_short=0.0d0
+ time_escp_long=0.0d0
+ time_escpsetup=0.0d0
+ time_escpcalc=0.0d0
#endif
cd print *,' in SET_TIMERS stime=',stime
return
include 'COMMON.IOUNITS'
include 'COMMON.TIME1'
include 'COMMON.SETUP'
+ include 'COMMON.MD'
#ifdef MPI
time1=MPI_WTIME()
write (iout,'(80(1h=)/a/(80(1h=)))')
& time_bcast+time_reduce+time_gather+time_scatter+
& time_sendrecv+time_barrier_g+time_barrier_e+time_bcastc
write (*,*) "Processor",fg_rank,myrank," enecalc",time_enecalc
+#ifdef TIMING_ENE
+ if (RESPA) then
+ write (*,*) "Processor",fg_rank,myrank," evdw_long",
+ & time_evdw_long
+ write (*,*) "Processor",fg_rank,myrank," evdw_short",
+ & time_evdw_short
+ write (*,*) "Processor",fg_rank,myrank," eelec_long",
+ & time_eelec_long
+ write (*,*) "Processor",fg_rank,myrank," eelec_short",
+ & time_eelec_short
+ write (*,*) "Processor",fg_rank,myrank," escp_long",
+ & time_escp_long
+ write (*,*) "Processor",fg_rank,myrank," escp_short",
+ & time_escp_short
+ else
+ write (*,*) "Processor",fg_rank,myrank," evdw",time_evdw
+ write (*,*) "Processor",fg_rank,myrank," eelec",time_eelec
+ write (*,*) "Processor",fg_rank,myrank," escp",time_escp
+ write (*,*) "Processor",fg_rank,myrank," escpsetup",
+ & time_escpsetup
+ write (*,*) "Processor",fg_rank,myrank," escpcalc",
+ & time_escpcalc
+ endif
+#endif
write (*,*) "Processor",fg_rank,myrank," sumene",time_sumene
write (*,*) "Processor",fg_rank,myrank," intfromcart",
& time_intfcart
include 'COMMON.REMD'
include 'COMMON.MD'
include 'COMMON.SBRIDGE'
- integer i,icall,iretcode,nfun
+ integer i,it,icall,iretcode,nfun
common /srutu/ icall
integer nharp,iharp(4,maxres/3)
integer nft_sc
- logical fail
+ logical fail,secondary_str /.true./
double precision energy(0:n_ene),etot,etota
double precision energy_long(0:n_ene),energy_short(0:n_ene)
double precision rms,frac,frac_nn,co
write (iout,*) "Energy evaluation/minimization"
call chainbuild_cart
c print *,'dc',dc(1,0),dc(2,0),dc(3,0)
+ if (nran_start.gt.0) then
+ write (iout,*)
+ & "Chains will be regenerated starting from residue",nran_start
+ do it=1,100
+ call gen_rand_conf_mchain(nran_start,*10)
+ write (iout,*) "Conformation successfully generated",it
+ goto 11
+ 10 write (iout,*) "Problems with regenerating chains",it
+ enddo
+ 11 continue
+ write (iout,*) "Cartesian coords after chain rebuild"
+ call cartprint
+ call chainbuild_cart
+ write (iout,*) "Cartesian coords after chainbuild_ecart"
+ call cartprint
+ call int_from_cart1(.false.)
+ call intout
+ endif
if (split_ene) then
print *,"Processor",myrank," after chainbuild"
icall=1
etota = energy(0)
etot =etota
call enerprint(energy(0))
+ if (secondary_str) then
call hairpin(.true.,nharp,iharp)
c print *,'after hairpin'
call secondary2(.true.)
+ endif
c print *,'after secondary'
if (minim) then
crc overlap test
write (iout,*) 'Calling OVERLAP_SC'
call overlap_sc(fail)
write (iout,*) "After overlap_sc"
+c cartname=prefix(:ilen(prefix))//'.x'
+c potE=etot
+c call cartoutx(0.0d0)
endif
if (searchsc) then
#endif
print *,'# eval/s',evals
print *,'refstr=',refstr
+ if (secondary_str) then
call hairpin(.false.,nharp,iharp)
- print *,'after hairpin'
+c print *,'after hairpin'
call secondary2(.true.)
- print *,'after secondary'
+c print *,'after secondary'
+ endif
call etotal(energy(0))
etot = energy(0)
call enerprint(energy(0))
projects / unres.git / commitdiff