--- /dev/null
+-1.0 5.0
+-1.0 5.0
+-1.0 5.0
+-1.0 5.0
+-1.0 5.0
itordp_nucl= 130
! ielep_nucl= 131
isidep_nucl=132
+ iscpp_nucl=133
!el local variables
integer :: ind_scint=0,ind_scint_old,ii,jj,i,j,iint,itmp
-
+ integer :: ind_scint_nucl=0
#ifdef MPI
integer :: my_sc_int(0:nfgtasks-1),my_ele_int(0:nfgtasks-1)
integer :: my_sc_intt(0:nfgtasks),my_ele_intt(0:nfgtasks)
integer :: n_sc_int_tot,my_sc_inde,my_sc_inds,ind_sctint,npept
+ integer :: n_sc_int_tot_nucl,my_sc_inde_nucl,my_sc_inds_nucl, &
+ ind_sctint_nucl,npept_nucl
+
integer :: nele_int_tot,my_ele_inds,my_ele_inde,ind_eleint_old,&
ind_eleint,ijunk,nele_int_tot_vdw,my_ele_inds_vdw,&
my_ele_inde_vdw,ind_eleint_vdw,ind_eleint_vdw_old,&
ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,&
ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,&
ichunk,int_index_old
+ integer :: nele_int_tot_nucl,my_ele_inds_nucl,my_ele_inde_nucl,&
+ ind_eleint_old_nucl,ind_eleint_nucl,nele_int_tot_vdw_nucl,&
+ my_ele_inds_vdw_nucl,my_ele_inde_vdw_nucl,ind_eleint_vdw_nucl,&
+ ind_eleint_vdw_old_nucl,nscp_int_tot_nucl,my_scp_inds_nucl,&
+ my_scp_inde_nucl,ind_scpint_nucl,ind_scpint_old_nucl
integer,dimension(5) :: nct_molec,nnt_molec
!el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
!el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
iatsc_e=nct-1
#endif
if (iatsc_s.eq.0) iatsc_s=1
+!----------------- scaling for nucleic acid GB
+ n_sc_int_tot_nucl=(nct_molec(2)-nnt_molec(2)+1)*(nct_molec(2)-nnt_molec(2))/2
+ call int_bounds(n_sc_int_tot_nucl,my_sc_inds_nucl,my_sc_inde_nucl)
+!write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
+ if (lprint) &
+ write (iout,*) 'Processor',fg_rank,' CG group',kolor,&
+ ' absolute rank',MyRank,&
+ ' n_sc_int_tot',n_sc_int_tot_nucl,' my_sc_inds=',my_sc_inds_nucl,&
+ ' my_sc_inde',my_sc_inde_nucl
+ ind_sctint_nucl=0
+ iatsc_s_nucl=0
+ iatsc_e_nucl=0
+ do i=1,nres !el !maxres
+ nint_gr_nucl(i)=0
+ nscp_gr_nucl(i)=0
+ ielstart_nucl(i)=0
+ ielend_nucl(i)=0
+ do j=1,maxint_gr
+ istart_nucl(i,j)=0
+ iend_nucl(i,j)=0
+ iscpstart_nucl(i,j)=0
+ iscpend_nucl(i,j)=0
+ enddo
+ enddo
+ do i=nnt_molec(2),nct_molec(2)-1
+! print*, "inloop2",i
+ call int_partition(ind_scint_nucl,my_sc_inds_nucl,my_sc_inde_nucl,i,&
+ iatsc_s_nucl,iatsc_e_nucl,i+1,nct_molec(2),nint_gr_nucl(i), &
+ istart_nucl(i,1),iend_nucl(i,1),*112)
+ print *,istart_nucl(i,1)
+ enddo
+ 112 continue
+ if (iatsc_s_nucl.eq.0) iatsc_s_nucl=1
+ print *,"tu mam",iatsc_s_nucl,iatsc_e_nucl
+
#ifdef MPI
if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,&
' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
enddo
endif
! lprint=.false.
+ write (iout,'(a)') 'Interaction array2:'
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ write (iout,'(i3,2(2x,2i4))') &
+ i,(istart_nucl(i,iint),iend_nucl(i,iint),iint=1,nint_gr_nucl(i))
+ enddo
+
ispp=4 !?? wham ispp=2
#ifdef MPI
! Now partition the electrostatic-interaction array
enddo ! i
13 continue
if (iatel_s.eq.0) iatel_s=1
+!----------now nucleic acid
+! if (itype(nres_molec(2),2).eq.ntyp1_molec(2)) then
+ npept_nucl=nct_molec(2)-nnt_molec(2)
+! else
+! npept_nucl=nct_molec(2)-nnt_molec(2)
+! endif
+ nele_int_tot_nucl=(npept_nucl-ispp)*(npept_nucl-ispp+1)/2
+ call int_bounds(nele_int_tot_nucl,my_ele_inds_nucl,my_ele_inde_nucl)
+ if (lprint) &
+ write (*,*) 'Processor',fg_rank,' CG group',kolor,&
+ ' absolute rank',MyRank,&
+ ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,&
+ ' my_ele_inde',my_ele_inde
+ iatel_s_nucl=0
+ iatel_e_nucl=0
+ ind_eleint_nucl=0
+ ind_eleint_old_nucl=0
+! if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+! nct_molec(1)=nres_molec(1)-1
+! else
+! nct_molec(1)=nres_molec(1)
+! endif
+! print *,"nct",nct,nct_molec(1),itype(nres_molec(1),1),ntyp_molec(1)
+ do i=nnt_molec(2),nct_molec(2)-3
+ ijunk=0
+ call int_partition(ind_eleint_nucl,my_ele_inds_nucl,my_ele_inde_nucl,i,&
+ iatel_s_nucl,iatel_e_nucl,i+ispp,nct_molec(2)-1,&
+ ijunk,ielstart_nucl(i),ielend_nucl(i),*113)
+ enddo ! i
+ 113 continue
+ if (iatel_s.eq.0) iatel_s=1
+
nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
! 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)
enddo ! i
if (iatel_s_vdw.eq.0) iatel_s_vdw=1
15 continue
+ if (iatel_s.eq.0) iatel_s=1
+
+ nele_int_tot_vdw_nucl=(npept_nucl-2)*(npept_nucl-2+1)/2
+! write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
+ call int_bounds(nele_int_tot_vdw_nucl,my_ele_inds_vdw_nucl,&
+ my_ele_inde_vdw_nucl)
+! write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
+! & " my_ele_inde_vdw",my_ele_inde_vdw
+ ind_eleint_vdw_nucl=0
+ ind_eleint_vdw_old_nucl=0
+ iatel_s_vdw_nucl=0
+ iatel_e_vdw_nucl=0
+ do i=nnt_molec(2),nct_molec(2)-3
+ ijunk=0
+ call int_partition(ind_eleint_vdw_nucl,my_ele_inds_vdw_nucl,&
+ my_ele_inde_vdw_nucl,i,&
+ iatel_s_vdw_nucl,iatel_e_vdw_nucl,i+2,nct_molec(2)-1,&
+ ijunk,ielstart_vdw_nucl(i),&
+ ielend_vdw(i),*115)
+! write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
+! & " ielend_vdw",ielend_vdw(i)
+ enddo ! i
+ if (iatel_s_vdw.eq.0) iatel_s_vdw=1
+ 115 continue
+
#else
iatel_s=nnt
iatel_e=nct_molec(1)-5 ! ?? wham iatel_e=nct-3
endif ! lprint
! iscp=3
iscp=2
+ iscp_nucl=2
! Partition the SC-p interaction array
#ifdef MPI
nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
endif
enddo ! i
14 continue
+ print *,"before inloop3",iatscp_s,iatscp_e,iscp_nucl
+ nscp_int_tot_nucl=(npept_nucl-iscp_nucl+1)*(npept_nucl-iscp_nucl+1)
+ call int_bounds(nscp_int_tot_nucl,my_scp_inds_nucl,my_scp_inde_nucl)
+ if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,&
+ ' absolute rank',myrank,&
+ ' nscp_int_tot',nscp_int_tot_nucl,' my_scp_inds=',my_scp_inds_nucl,&
+ ' my_scp_inde',my_scp_inde_nucl
+ print *,"nscp_int_tot_nucl",nscp_int_tot_nucl,my_scp_inds_nucl,my_scp_inde_nucl
+ iatscp_s_nucl=0
+ iatscp_e_nucl=0
+ ind_scpint_nucl=0
+ ind_scpint_old_nucl=0
+ do i=nnt_molec(2),nct_molec(2)-1
+ print *,"inloop3",i,nnt_molec(2)+iscp,nct_molec(2)-iscp
+ if (i.lt.nnt_molec(2)+iscp) then
+!d write (iout,*) 'i.le.nnt+iscp'
+ call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+ my_scp_inde_nucl,i,iatscp_s_nucl,iatscp_e_nucl,i+iscp,&
+ nct_molec(2),nscp_gr_nucl(i),iscpstart_nucl(i,1),&
+ iscpend_nucl(i,1),*114)
+ else if (i.gt.nct_molec(2)-iscp) then
+!d write (iout,*) 'i.gt.nct-iscp'
+ call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+ my_scp_inde_nucl,i,&
+ iatscp_s_nucl,iatscp_e_nucl,nnt_molec(2),i-iscp,nscp_gr_nucl(i),&
+ iscpstart_nucl(i,1),&
+ iscpend_nucl(i,1),*114)
+ else
+ call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+ my_scp_inde_nucl,i,iatscp_s_nucl,iatscp_e_nucl,nnt_molec(2),&
+ i-iscp,nscp_gr_nucl(i),iscpstart_nucl(i,1),&
+ iscpend_nucl(i,1),*114)
+ ii=nscp_gr_nucl(i)+1
+ call int_partition(ind_scpint_nucl,my_scp_inds_nucl,&
+ my_scp_inde_nucl,i,iatscp_s_nucl,iatscp_e_nucl,i+iscp,&
+ nct_molec(2),nscp_gr_nucl(i),iscpstart_nucl(i,ii),&
+ iscpend_nucl(i,ii),*114)
+ endif
+ enddo ! i
+ 114 continue
+ print *, "after inloop3",iatscp_s_nucl,iatscp_e_nucl
#else
iatscp_s=nnt
iatscp_e=nct_molec(1)-1
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
+! do i=igrad_start,igrad_end
+! write(*,*) 'Processor:',fg_rank,myrank,i,&
+! jgrad_start(i),jgrad_end(i)
+! enddo
endif
if (nfgtasks.gt.1) then
call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,&
real(kind=8),dimension(3) :: erij,gg,gg_lipi,gg_lipj
!-----------------------------------------------------------------------------
end module calc_data
+ module calc_data_nucl
+!-----------------------------------------------------------------------------
+! commom.calc common/calc/
+ integer :: i,j,k,l,ind
+ real(kind=8) :: rrij,rij,xj,yj,zj,dxj,dyj,dzj,&
+ chi1,chi2,chi12,chip1,chip2,chip12,alf1,alf2,alf12,om1,om2,om12,&
+ om1om2,chiom1,chiom2,chiom12,chipom1,chipom2,chipom12,eps1,&
+ faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,&
+ sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,&
+ eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,&
+ sigder,dsci_inv,dscj_inv
+ real(kind=8),dimension(3) :: erij,gg,gg_lipi,gg_lipj
+!-----------------------------------------------------------------------------
+ end module calc_data_nucl
+
#endif
real(kind=8),dimension(:),allocatable :: weights !(n_ene)
real(kind=8) :: temp0,scal14,cutoff_corr,delt_corr,r0_corr
- integer :: ipot
+ integer :: ipot,ipot_nucl
! common /potentials/
character(len=3),dimension(5) :: potname = &
(/'LJ ','LJK','BP ','GB ','GBV'/)
real(kind=8),dimension(:,:),allocatable :: aa_aq,bb_aq,augm,aa_lip,bb_lip !(ntyp,ntyp)
real(kind=8),dimension(:),allocatable :: sc_aa_tube_par,sc_bb_tube_par,&
acavtub,bcavtub,ccavtub,dcavtub,tubetranene
+ real(kind=8),dimension(:,:),allocatable :: aa_nucl,bb_nucl
real(kind=8) :: acavtubpep,bcavtubpep,ccavtubpep,dcavtubpep, &
tubetranenepep,pep_aa_tube,pep_bb_tube,tubeR0
real(kind=8),dimension(3) :: tubecenter
real(kind=8),dimension(:,:),allocatable :: aad,bad !(ntyp,2)
real(kind=8),dimension(2,2) :: app,bpp,ael6,ael3
+ real(kind=8),dimension(:),allocatable :: aad_nucl,bad_nucl !(ntyp,2)
+ real(kind=8),dimension(2,2) :: app_nucl,bpp_nucl
+ real(kind=8),dimension(:,:),allocatable :: ael6_nucl,&
+ ael3_nucl,ael32_nucl,ael63_nucl
integer :: expon,expon2, nnt,nct,itypro
integer,dimension(:,:),allocatable :: istart,iend !(maxres,maxint_gr)
integer,dimension(:),allocatable :: nint_gr,itel,&
ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr !(maxres)
+ integer,dimension(:,:),allocatable :: istart_nucl,iend_nucl !(maxres,maxint_gr)
+ integer,dimension(:),allocatable :: nint_gr_nucl,itel_nucl,&
+ ielstart_nucl,ielend_nucl,ielstart_vdw_nucl,ielend_vdw_nucl,nscp_gr_nucl !(maxres)
+ integer,dimension(:,:),allocatable :: iscpstart_nucl,iscpend_nucl !(maxres,maxint_gr)
+
integer,dimension(:),allocatable :: istype,molnum
integer,dimension(:,:),allocatable :: itype ! now itype has more molecule types
integer,dimension(:,:),allocatable :: iscpstart,iscpend !(maxres,maxint_gr)
integer :: iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,&
iatel_e_vdw,iatscp_s,iatscp_e,ispp,iscp
+ integer :: iatsc_s_nucl,iatsc_e_nucl,iatel_s_nucl,iatel_e_nucl,&
+ iatel_s_vdw_nucl,iatel_e_vdw_nucl,iatscp_s_nucl,iatscp_e_nucl,&
+ ispp_nucl,iscp_nucl
+
! 12/1/95 Array EPS included in the COMMON block.
! common /body/
real(kind=8),dimension(:,:),allocatable :: sigma !(0:ntyp1,0:ntyp1)
real(kind=8),dimension(:),allocatable :: chip,alp,sigma0,&
sigii,rr0 !(ntyp)
real(kind=8),dimension(2,2) :: rpp,epp,elpp6,elpp3
+ real(kind=8),dimension(:,:),allocatable :: sigma_nucl !(0:ntyp1,0:ntyp1)
+ real(kind=8),dimension(:,:),allocatable :: eps_nucl,sigmaii_nucl,&
+ chi_nucl,r0_nucl, chip_nucl !(ntyp,ntyp) r0e !!! nie używane
+ real(kind=8),dimension(:),allocatable :: alp_nucl,sigma0_nucl,&
+ sigii_nucl,rr0_nucl !(ntyp)
+ real(kind=8),dimension(2,2) :: rpp_nucl,epp_nucl
+ real(kind=8),dimension(:,:),allocatable ::elpp6_nucl,&
+ elpp3_nucl,elpp32_nucl,elpp63_nucl
+ real(kind=8):: r0pp,epspp,AEES,BEES
+
real(kind=8),dimension(:,:),allocatable :: r0d,eps_scp,rscp !(ntyp,2) r0d !!! nie używane
+ real(kind=8),dimension(:),allocatable :: eps_scp_nucl,rscp_nucl!(ntyp,2) r0d !!! nie używane
+
! 12/5/03 modified 09/18/03 Bond stretching parameters.
! common /stretch/
real(kind=8) :: vbldp0,akp,distchainmax,vbldpDUM
ientout,izs1,isecpred,ibond,irest2,iifrag,icart,irest1,isccor,&
ithep_pdb,irotam_pdb,iliptranpar,itube, &
ibond_nucl,ithep_nucl,irotam_nucl,itorp_nucl,itordp_nucl, &
- isidep_nucl
+ isidep_nucl,iscpp_nucl
#ifdef WHAM_RUN
! el wham iounits
integer :: isidep1,ihist,iweight,izsc,idistr
fouriername,elename,sidename,scpname,sccorname,patname,&
thetname_pdb,rotname_pdb,liptranname,tubename, &
bondname_nucl,thetname_nucl,rotname_nucl,torname_nucl, &
- tordname_nucl,sidename_nucl
+ tordname_nucl,sidename_nucl,scpname_nucl
!-----------------------------------------------------------------------
! INP - main input file
! IOUT - list file
(/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25,&
26,27,28,29,30,31,32,33,34,35,36,37,38/)
-
+ character(len=1), dimension(2) :: sugartyp = (/'D',' '/)
!#endif
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,&
grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres)
!-----------------------------NUCLEIC GRADIENT
- real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl
+ real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl, &
+ gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc
! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
call ebond_nucl(estr_nucl)
call ebend_nucl(ebe_nucl)
call etor_nucl(etors_nucl)
+ call esb_gb(evdwsb,eelsb)
+! call multibody_hb(ecorr,ecorr3,n_corr,n_corr1)
+ call epp_nucl_sub(evdwpp,eespp)
+ call epsb(evdwpsb,eelpsb)
+
print *,"after ebend", ebe_nucl
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
iteli=itel(i)
itelj=itel(j)
if (j.eq.i+2 .and. itelj.eq.2) iteli=2
- aaa=app(iteli,itelj)
- bbb=bpp(iteli,itelj)
- ael6i=ael6(iteli,itelj)
- ael3i=ael3(iteli,itelj)
+ aaa=app_nucl(iteli,itelj)
+ bbb=bpp_nucl(iteli,itelj)
+ ael6i=ael6_nucl(iteli,itelj)
+ ael3i=ael3_nucl(iteli,itelj)
dxj=dc(1,j)
dyj=dc(2,j)
dzj=dc(3,j)
allocate(ielstart_vdw(nres))
allocate(ielend_vdw(nres))
!(maxres)
+ allocate(nint_gr_nucl(nres))
+ allocate(nscp_gr_nucl(nres))
+ allocate(ielstart_nucl(nres))
+ allocate(ielend_nucl(nres))
+!(maxres)
+ allocate(istart_nucl(nres,maxint_gr))
+ allocate(iend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(iscpstart_nucl(nres,maxint_gr))
+ allocate(iscpend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(ielstart_vdw_nucl(nres))
+ allocate(ielend_vdw_nucl(nres))
allocate(lentyp(0:nfgtasks-1))
!(0:maxprocs-1)
allocate(gradafm(3,-1:nres))
allocate(gradb_nucl(3,-1:nres))
allocate(gradbx_nucl(3,-1:nres))
+ allocate(gvdwpsb1(3,-1:nres))
+ allocate(gelpp(3,-1:nres))
+ allocate(gvdwpsb(3,-1:nres))
+ allocate(gelsbc(3,-1:nres))
+ allocate(gelsbx(3,-1:nres))
+ allocate(gvdwsbx(3,-1:nres))
+ allocate(gvdwsbc(3,-1:nres))
!(3,maxres)
allocate(grad_shield_side(3,50,nres))
allocate(grad_shield_loc(3,50,nres))
real(kind=8) :: difi,thetiii
integer itheta
etheta_nucl=0.0D0
- print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres
+! print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres
do i=ithet_nucl_start,ithet_nucl_end
if ((itype(i-1,2).eq.ntyp1_molec(2)).or.&
(itype(i-2,2).eq.ntyp1_molec(2)).or. &
i,theta(i)*rad2deg,phii*rad2deg, &
phii1*rad2deg,ethetai
etheta_nucl=etheta_nucl+ethetai
- print *,i,"partial sum",etheta_nucl
+! print *,i,"partial sum",etheta_nucl
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1
gloc(nphi+i-2,icg)=wang_nucl*dethetai
lprn=.false.
! lprn=.true.
etors_nucl=0.0D0
- print *,"iphi_nucl_start/end", iphi_nucl_start,iphi_nucl_end
+! print *,"iphi_nucl_start/end", iphi_nucl_start,iphi_nucl_end
do i=iphi_nucl_start,iphi_nucl_end
if (itype(i-2,2).eq.ntyp1_molec(2) .or. itype(i-1,2).eq.ntyp1_molec(2) &
.or. itype(i-3,2).eq.ntyp1_molec(2) &
itori=itortyp_nucl(itype(i-2,2))
itori1=itortyp_nucl(itype(i-1,2))
phii=phi(i)
- print *,i,itori,itori1
+! print *,i,itori,itori1
gloci=0.0D0
!C Regular cosine and sine terms
do j=1,nterm_nucl(itori,itori1)
enddo
return
end subroutine etor_nucl
+!------------------------------------------------------------
+ subroutine epp_nucl_sub(evdw1,ees)
+!C
+!C This subroutine calculates the average interaction energy and its gradient
+!C in the virtual-bond vectors between non-adjacent peptide groups, based on
+!C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+!C The potential depends both on the distance of peptide-group centers and on
+!C the orientation of the CA-CA virtual bonds.
+!C
+ integer :: i,j,k,iteli,itelj,num_conti,isubchap,ind
+ real(kind=8) :: dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
+ real(kind=8) :: xj,yj,zj,rij,rrmij,sss,r3ij,r6ij,evdw1,&
+ dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
+ dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,sss_grad,fac,evdw1ij
+ integer xshift,yshift,zshift
+ real(kind=8),dimension(3):: ggg,gggp,gggm,erij
+ real(kind=8) :: ees,eesij
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ ind=0
+!c
+!c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
+!c
+ print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+ do i=iatel_s_nucl,iatel_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) 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
+
+ do j=ielstart_nucl(i),ielend_nucl(i)
+ if (itype(j,2).eq.ntyp1_molec(2) .or. itype(j+1,2).eq.ntyp1_molec(2)) cycle
+ ind=ind+1
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(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
+ isubchap=0
+ 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
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+!C print *,i,j
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
+
+ rij=xj*xj+yj*yj+zj*zj
+!c write (2,*)"ij",i,j," r0pp",r0pp," rij",rij," epspp",epspp
+ fac=(r0pp**2/rij)**3
+ ev1=epspp*fac*fac
+ ev2=epspp*fac
+ evdw1ij=ev1-2*ev2
+ fac=(-ev1-evdw1ij)/rij
+! write (2,*)"fac",fac," ev1",ev1," ev2",ev2," evdw1ij",evdw1ij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"evdw1ij",evdw1ij
+ evdw1=evdw1+evdw1ij
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+ do k=1,3
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ enddo
+!c phoshate-phosphate electrostatic interactions
+ rij=dsqrt(rij)
+ fac=1.0d0/rij
+ eesij=dexp(-BEES*rij)*fac
+! write (2,*)"fac",fac," eesijpp",eesij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"eesijpp",eesij
+ ees=ees+eesij
+!c fac=-eesij*fac
+ fac=-(fac+BEES)*eesij*fac
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+!c write(2,*) "ggg",i,j,ggg(1),ggg(2),ggg(3)
+!c write(2,*) "gelpp",i,(gelpp(k,i),k=1,3)
+!c write(2,*) "gelpp",j,(gelpp(k,j),k=1,3)
+ do k=1,3
+ gelpp(k,i)=gelpp(k,i)-ggg(k)
+ gelpp(k,j)=gelpp(k,j)+ggg(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!c ees=332.0d0*ees
+ ees=AEES*ees
+ do i=nnt,nct
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ do k=1,3
+ gvdwpp(k,i)=6*gvdwpp(k,i)
+!c gelpp(k,i)=332.0d0*gelpp(k,i)
+ gelpp(k,i)=AEES*gelpp(k,i)
+ enddo
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ enddo
+!c write (2,*) "total EES",ees
+ return
+ end subroutine epp_nucl_sub
+!---------------------------------------------------------------------
+ subroutine epsb(evdwpsb,eelpsb)
+! use comm_locel
+!C
+!C This subroutine calculates the excluded-volume interaction energy between
+!C peptide-group centers and side chains and its gradient in virtual-bond and
+!C side-chain vectors.
+!C
+ real(kind=8),dimension(3):: ggg
+ integer :: i,iint,j,k,iteli,itypj,subchap
+ real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
+ e1,e2,evdwij,rij,evdwpsb,eelpsb
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
+
+!cd print '(a)','Enter ESCP'
+!cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+ eelpsb=0.0d0
+ evdwpsb=0.0d0
+ print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl
+ do i=iatscp_s_nucl,iatscp_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) &
+ .or. itype(i+1,2).eq.ntyp1_molec(2)) 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))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do iint=1,nscp_gr_nucl(i)
+
+ do j=iscpstart_nucl(i,iint),iscpend_nucl(i,iint)
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+!C Uncomment following three lines for SC-p interactions
+!c xj=c(1,nres+j)-xi
+!c yj=c(2,nres+j)-yi
+!c zj=c(3,nres+j)-zi
+!C 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)
+ 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-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ 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-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ fac=rrij**expon2
+ e1=fac*fac*aad_nucl(itypj)
+ e2=fac*bad_nucl(itypj)
+ if (iabs(j-i) .le. 2) then
+ e1=scal14*e1
+ e2=scal14*e2
+ endif
+ evdwij=e1+e2
+ evdwpsb=evdwpsb+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a4)') &
+ 'evdw2',i,j,evdwij,"tu4"
+!C
+!C Calculate contributions to the gradient in the virtual-bond and SC vectors.
+!C
+ fac=-(evdwij+e1)*rrij
+ ggg(1)=xj*fac
+ ggg(2)=yj*fac
+ ggg(3)=zj*fac
+ do k=1,3
+ gvdwpsb1(k,i)=gvdwpsb1(k,i)-ggg(k)
+ gvdwpsb(k,j)=gvdwpsb(k,j)+ggg(k)
+ enddo
+ enddo
+
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwpsb(j,i)=expon*gvdwpsb(j,i)
+ gvdwpsb1(j,i)=expon*gvdwpsb1(j,i)
+ enddo
+ enddo
+ return
+ end subroutine epsb
+
+!------------------------------------------------------
+ subroutine esb_gb(evdwsb,eelsb)
+ use comm_locel
+ use calc_data_nucl
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0iji,evdwsb,eelsb,ecorr,eelij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,faclip,sig0ij
+ integer :: ii
+ logical lprn
+ evdw=0.0D0
+ eelsb=0.0d0
+ ecorr=0.0d0
+ evdwsb=0.0D0
+ lprn=.false.
+ ind=0
+! print *,"iastsc_nucl",iatsc_s_nucl,iatsc_e_nucl
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ num_conti=0
+ itypi=itype(i,2)
+! PRINT *,"I=",i,itypi
+ if (itypi.eq.ntyp1_molec(2)) cycle
+ itypi1=itype(i+1,2)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+!C
+!C Calculate SC interaction energy.
+!C
+ do iint=1,nint_gr_nucl(i)
+! print *,"tu?",i,istart_nucl(i,iint),iend_nucl(i,iint)
+ do j=istart_nucl(i,iint),iend_nucl(i,iint)
+ ind=ind+1
+! print *,"JESTEM"
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma_nucl(itypi,itypj)
+ chi1=chi_nucl(itypi,itypj)
+ chi2=chi_nucl(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip_nucl(itypi,itypj)
+ chip2=chip_nucl(itypj,itypi)
+ chip12=chip1*chip2
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ 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-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+!C Calculate angle-dependent terms of energy and contributions to their
+!C derivatives.
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular_nucl
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+! print *,rij_shift,"rij_shift"
+!c write (2,*) " rij",1.0D0/rij," sig",sig," sig0ij",sig0ij,
+!c & " rij_shift",rij_shift
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+!c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa_nucl(itypi,itypj)
+ e2=fac*bb_nucl(itypi,itypj)
+ evdwij=eps1*eps2rt*(e1+e2)
+!c write (2,*) "eps1",eps1," eps2rt",eps2rt,
+!c & " e1",e1," e2",e2," evdwij",evdwij
+ eps2der=evdwij
+ evdwij=evdwij*eps2rt
+ evdwsb=evdwsb+evdwij
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,2),i,restyp(itypj,2),j, &
+ epsi,sigm,chi1,chi2,chip1,chip2, &
+ eps1,eps2rt**2,sig,sig0ij, &
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij
+ write (iout,*) "aa",aa_nucl(itypi,itypj)," bb",bb_nucl(itypi,itypj)
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,e15.3,a4)') &
+ 'evdw',i,j,evdwij,"tu3"
+
+
+!C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+!c fac=0.0d0
+!C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+!C Calculate angular part of the gradient.
+ call sc_grad_nucl
+ call eelsbij(eelij)
+ if (energy_dec .and. &
+ (j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
+ write (istat,'(e14.5)') evdwij
+ eelsb=eelsb+eelij
+ enddo ! j
+ enddo ! iint
+ num_cont_hb(i)=num_conti
+ enddo ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!cccc energy_dec=.false.
+ return
+ end subroutine esb_gb
+!-------------------------------------------------------------------------------
+ subroutine eelsbij(eesij)
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,dcosb,dcosg
+ real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,rlocshield,fracinbuf
+ integer xshift,yshift,zshift,ilist,iresshield
+
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ integer :: iteli,itelj,kkk,kkll,m,isubchap
+ real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp,facfac
+ real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i,ael63i,ael32i
+ real(kind=8) :: dx_normj,dy_normj,dz_normj,&
+ r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,fac5,fac6,&
+ el1,el2,el3,el4,eesij,ees0ij,facvdw,facel,fac1,ecosa,&
+ ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,&
+ a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,&
+ ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,&
+ ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,&
+ ecosgp,ecosam,ecosbm,ecosgm,ghalf,itypi,itypj
+ ind=ind+1
+ itypi=itype(i,2)
+ itypj=itype(j,2)
+! print *,i,j,itypi,itypj,istype(i),istype(j),"????"
+ ael6i=ael6_nucl(itypi,itypj)
+ ael3i=ael3_nucl(itypi,itypj)
+ ael63i=ael63_nucl(itypi,itypj)
+ ael32i=ael32_nucl(itypi,itypj)
+!c write (iout,*) "eelecij",i,j,itype(i),itype(j),
+!c & ael6i,ael3i,ael63i,al32i,rij,rrij
+ dxj=dc(1,j+nres)
+ dyj=dc(2,j+nres)
+ dzj=dc(3,j+nres)
+ dx_normi=dc_norm(1,i+nres)
+ dy_normi=dc_norm(2,i+nres)
+ dz_normi=dc_norm(3,i+nres)
+ dx_normj=dc_norm(1,j+nres)
+ dy_normj=dc_norm(2,j+nres)
+ dz_normj=dc_norm(3,j+nres)
+!c xj=c(1,j)+0.5D0*dxj-xmedi
+!c yj=c(2,j)+0.5D0*dyj-ymedi
+!c zj=c(3,j)+0.5D0*dzj-zmedi
+ if (ipot_nucl.ne.2) then
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ else
+ cosa=om12
+ cosb=om1
+ cosg=om2
+ endif
+ r3ij=rij*rrij
+ r6ij=r3ij*r3ij
+ fac=cosa-3.0D0*cosb*cosg
+ facfac=fac*fac
+ fac1=3.0d0*(cosb*cosb+cosg*cosg)
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ fac5=ael63i*r6ij
+ fac6=ael32i*r6ij
+!c write (iout,*) "r3ij",r3ij," r6ij",r6ij," fac",fac," fac1",fac1,
+!c & " fac2",fac2," fac3",fac3," fac4",fac4," fac5",fac5," fac6",fac6
+ el1=fac3*(4.0D0+facfac-fac1)
+ el2=fac4*fac
+ el3=fac5*(2.0d0-2.0d0*facfac+fac1)
+ el4=fac6*facfac
+ eesij=el1+el2+el3+el4
+!C 12/26/95 - for the evaluation of multi-body H-bonding interactions
+ ees0ij=4.0D0+facfac-fac1
+
+ if (energy_dec) then
+ if(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2) &
+ write (istat,'(2a1,i4,1x,2a1,i4,4f10.5,3e12.5,$)') &
+ sugartyp(istype(i)),restyp(itypi,2),i,sugartyp(istype(j)),&
+ restyp(itypj,2),j,1.0d0/rij,cosa,cosb,cosg,fac*r3ij, &
+ (4.0D0+facfac-fac1)*r6ij,(2.0d0-2.0d0*facfac+fac1)*r6ij
+ write (iout,'(a6,2i5,e15.3)') 'ees',i,j,eesij
+ endif
+
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ facel=-3.0d0*rrij*(eesij+el1+el3+el4)
+ fac1=fac
+!c erij(1)=xj*rmij
+!c erij(2)=yj*rmij
+!c erij(3)=zj*rmij
+!*
+!* Radial derivatives. First process both termini of the fragment (i,j)
+!*
+ ggg(1)=facel*xj
+ ggg(2)=facel*yj
+ ggg(3)=facel*zj
+ do k=1,3
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k)
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k)
+ enddo
+!*
+!* Angular part
+!*
+ ecosa=2.0D0*fac3*fac1+fac4+(-4.0d0*fac5+2.0d0*fac6)*fac1
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ fac5= 6.0d0*fac5
+ fac6=-6.0d0*fac6
+ ecosb=fac3*(fac1*cosg+cosb)+cosg*fac4+(cosb+2*fac1*cosg)*fac5+&
+ fac6*fac1*cosg
+ ecosg=fac3*(fac1*cosb+cosg)+cosb*fac4+(cosg+2*fac1*cosb)*fac5+&
+ fac6*fac1*cosb
+ do k=1,3
+ dcosb(k)=rij*(dc_norm(k,i+nres)-erij(k)*cosb)
+ dcosg(k)=rij*(dc_norm(k,j+nres)-erij(k)*cosg)
+ enddo
+ do k=1,3
+ ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ enddo
+ do k=1,3
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k) &
+ +(ecosa*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres))&
+ + ecosb*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k) &
+ +(ecosa*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosg*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ enddo
+ IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and. j.gt.i+1 .and.&
+ num_conti.le.maxconts) THEN
+!C
+!C Calculate the contact function. The ith column of the array JCONT will
+!C contain the numbers of atoms that make contacts with the atom I (of numbers
+!C greater than I). The arrays FACONT and GACONT will contain the values of
+!C the contact function and its derivative.
+ r0ij=2.20D0*sigma(itypi,itypj)
+!c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij
+ call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont)
+!c write (2,*) "fcont",fcont
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ if (num_conti.gt.maxconts) then
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',&
+ ' will skip next contacts for this conf.'
+ else
+ jcont_hb(num_conti,i)=j
+!c write (iout,*) "num_conti",num_conti,
+!c & " jcont_hb",jcont_hb(num_conti,i)
+!C Calculate contact energies
+ cosa4=4.0D0*cosa
+ wij=cosa-3.0D0*cosb*cosg
+ cosbg1=cosb+cosg
+ cosbg2=cosb-cosg
+ fac3=dsqrt(-ael6i)*r3ij
+!c write (2,*) "ael6i",ael6i," r3ij",r3ij," fac3",fac3
+ ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
+ if (ees0tmp.gt.0) then
+ ees0pij=dsqrt(ees0tmp)
+ else
+ ees0pij=0
+ endif
+ ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
+ if (ees0tmp.gt.0) then
+ ees0mij=dsqrt(ees0tmp)
+ else
+ ees0mij=0
+ endif
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+!c write (iout,*) "i",i," j",j,
+!c & " ees0m",ees0m(num_conti,i)," ees0p",ees0p(num_conti,i)
+ ees0pij1=fac3/ees0pij
+ ees0mij1=fac3/ees0mij
+ fac3p=-3.0D0*fac3*rrij
+ ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
+ ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
+ ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
+ ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
+ ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
+ ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
+ ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
+ ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
+ ecosap=ecosa1+ecosa2
+ ecosbp=ecosb1+ecosb2
+ ecosgp=ecosg1+ecosg2
+ ecosam=ecosa1-ecosa2
+ ecosbm=ecosb1-ecosb2
+ ecosgm=ecosg1-ecosg2
+!C End diagnostics
+ facont_hb(num_conti,i)=fcont
+ fprimcont=fprimcont/rij
+ do k=1,3
+ gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
+ gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
+ enddo
+ gggp(1)=gggp(1)+ees0pijp*xj
+ gggp(2)=gggp(2)+ees0pijp*yj
+ gggp(3)=gggp(3)+ees0pijp*zj
+ gggm(1)=gggm(1)+ees0mijp*xj
+ gggm(2)=gggm(2)+ees0mijp*yj
+ gggm(3)=gggm(3)+ees0mijp*zj
+!C Derivatives due to the contact function
+ gacont_hbr(1,num_conti,i)=fprimcont*xj
+ gacont_hbr(2,num_conti,i)=fprimcont*yj
+ gacont_hbr(3,num_conti,i)=fprimcont*zj
+ do k=1,3
+!c
+!c Gradient of the correlation terms
+!c
+ gacontp_hb1(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontp_hb2(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres)) &
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontp_hb3(k,num_conti,i)=gggp(k)
+ gacontm_hb1(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontm_hb2(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontm_hb3(k,num_conti,i)=gggm(k)
+ enddo
+ endif
+ endif
+ ENDIF
+ return
+ end subroutine eelsbij
+!------------------------------------------------------------------
+ subroutine sc_grad_nucl
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+sigder*sigsq_om2
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12+sigder*sigsq_om12
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwsbx(k,i)=gvdwsbx(k,i)-gg(k) &
+ +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwsbx(k,j)=gvdwsbx(k,j)+gg(k) &
+ +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!C
+!C Calculate the components of the gradient in DC and X
+!C
+ do l=1,3
+ gvdwsbc(l,i)=gvdwsbc(l,i)-gg(l)
+ gvdwsbc(l,j)=gvdwsbc(l,j)+gg(l)
+ enddo
+ return
+ end subroutine sc_grad_nucl
+
+!----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
end module energy
end subroutine sc_angular
!-----------------------------------------------------------------------------
! initialize_p.F
+ subroutine sc_angular_nucl
+! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
+! om12. Called by ebp, egb, and egbv.
+! use calc_data
+! implicit none
+! include 'COMMON.CALC'
+! include 'COMMON.IOUNITS'
+ use comm_locel
+ use calc_data_nucl
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ chiom12=chi12*om12
+! Calculate eps1(om12) and its derivative in om12
+ faceps1=1.0D0-om12*chiom12
+ faceps1_inv=1.0D0/faceps1
+ eps1=dsqrt(faceps1_inv)
+! Following variable is eps1*deps1/dom12
+ eps1_om12=faceps1_inv*chiom12
+! diagnostics only
+! faceps1_inv=om12
+! eps1=om12
+! eps1_om12=1.0d0
+! write (iout,*) "om12",om12," eps1",eps1
+! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
+! and om12.
+ om1om2=om1*om2
+ chiom1=chi1*om1
+ chiom2=chi2*om2
+ facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
+ sigsq=1.0D0-facsig*faceps1_inv
+ sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
+ sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
+ sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
+ chipom1=chip1*om1
+ chipom2=chip2*om2
+ chipom12=chip12*om12
+ facp=1.0D0-om12*chipom12
+ facp_inv=1.0D0/facp
+ facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
+! write (iout,*) "chipom1",chipom1," chipom2",chipom2,
+! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv
+! Following variable is the square root of eps2
+ eps2rt=1.0D0-facp1*facp_inv
+! Following three variables are the derivatives of the square root of eps
+! in om1, om2, and om12.
+ eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
+ eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
+ eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
+! Evaluate the "asymmetric" factor in the VDW constant, eps3
+ eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
+! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
+! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
+! & " eps2rt_om12",eps2rt_om12
+! Calculate whole angle-dependent part of epsilon and contributions
+! to its derivatives
+ return
+ end subroutine sc_angular_nucl
+
!-----------------------------------------------------------------------------
subroutine int_bounds(total_ints,lower_bound,upper_bound)
! implicit real*8 (a-h,o-z)
sccmj=0.0D0
do i=1,nscat
sccmj=sccmj+sccor(j,i)
- print *,"insccent", ires,sccor(j,i)
+!C print *,"insccent", ires,sccor(j,i)
enddo
dc(j,ires)=sccmj/nscat
enddo
endif
#endif
nct=nres
-!d print *,'NNT=',NNT,' NCT=',NCT
+ print *,'NNT=',NNT,' NCT=',NCT
if (itype(1,1).eq.ntyp1) nnt=2
if (itype(nres,1).eq.ntyp1) nct=nct-1
if (pdbref) then
#endif
allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
- print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
+! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
!el from energy module---------
allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
!el--------------------
read (itorp_nucl,*,end=113,err=113) &
(itortyp_nucl(i),i=1,ntyp_molec(2))
- print *,itortyp_nucl(:)
+! print *,itortyp_nucl(:)
!c write (iout,*) 'ntortyp',ntortyp
do i=1,ntortyp_nucl
do j=1,ntortyp_nucl
read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
- print *,nterm_nucl(i,j),nlor_nucl(i,j)
+! print *,nterm_nucl(i,j),nlor_nucl(i,j)
v0ij=0.0d0
si=-1.0d0
do k=1,nterm_nucl(i,j)
allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
+
allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
allocate(epslip(ntyp,ntyp))
end select
continue
close (isidep)
+
!-----------------------------------------------------------------------
! Calculate the "working" parameters of SC interactions.
endif
enddo
enddo
+
+ allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
+ allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
+
+! augm(:,:)=0.0D0
+! chip(:)=0.0D0
+! alp(:)=0.0D0
+! sigma0(:)=0.0D0
+! sigii(:)=0.0D0
+! rr0(:)=0.0D0
+
+ read (isidep_nucl,*) ipot_nucl
+! print *,"TU?!",ipot_nucl
+ if (ipot_nucl.eq.1) then
+ do i=1,ntyp_molec(2)
+ do j=i,ntyp_molec(2)
+ read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
+ elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
+ enddo
+ enddo
+ else
+ do i=1,ntyp_molec(2)
+ do j=i,ntyp_molec(2)
+ read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
+ chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
+ elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
+ enddo
+ enddo
+ endif
+! rpp(1,1)=2**(1.0/6.0)*5.16158
+ do i=1,ntyp_molec(2)
+ do j=i,ntyp_molec(2)
+ rrij=sigma_nucl(i,j)
+ r0_nucl(i,j)=rrij
+ r0_nucl(j,i)=rrij
+ rrij=rrij**expon
+ epsij=4*eps_nucl(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_nucl(i,j)=epsij*rrij*rrij
+ bb_nucl(i,j)=-sigeps*epsij*rrij
+ ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
+ ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
+ ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
+ ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
+ sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
+ 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
+ enddo
+ do j=1,i-1
+ aa_nucl(i,j)=aa_nucl(j,i)
+ bb_nucl(i,j)=bb_nucl(j,i)
+ ael3_nucl(i,j)=ael3_nucl(j,i)
+ ael6_nucl(i,j)=ael6_nucl(j,i)
+ ael63_nucl(i,j)=ael63_nucl(j,i)
+ ael32_nucl(i,j)=ael32_nucl(j,i)
+ elpp3_nucl(i,j)=elpp3_nucl(j,i)
+ elpp6_nucl(i,j)=elpp6_nucl(j,i)
+ elpp63_nucl(i,j)=elpp63_nucl(j,i)
+ elpp32_nucl(i,j)=elpp32_nucl(j,i)
+ eps_nucl(i,j)=eps_nucl(j,i)
+ sigma_nucl(i,j)=sigma_nucl(j,i)
+ sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
+ enddo
+ enddo
+
write(iout,*) "tube param"
read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
ccavtubpep,dcavtubpep,tubetranenepep
endif
! lprint=.false.
#endif
+ allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
+
+ do i=1,ntyp_molec(2)
+ read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
+ enddo
+ do i=1,ntyp_molec(2)
+ aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
+ bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
+ enddo
+ r0pp=1.12246204830937298142*5.16158
+ epspp=4.95713/4
+ AEES=108.661
+ BEES=0.433246
+
!
! Define the constants of the disulfide bridge
!
! Read the PDB file and convert the peptide geometry into virtual-chain
! geometry.
use geometry_data
- use energy_data, only: itype
+ use energy_data, only: itype,istype
use control_data
use compare_data
use MPI_data
real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
real(kind=8),dimension(:,:), allocatable :: c_temporary
integer,dimension(:,:) , allocatable :: itype_temporary
+ integer,dimension(:),allocatable :: istype_temp
efree_temp=0.0d0
ibeg=1
ishift1=0
ishift1=ishift1-1 !!!!!
! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
ires=ires-ishift+ishift1
- print *,ires,ishift,ishift1
+! print *,ires,ishift,ishift1
ires_old=ires
ibeg=0
else
molecule=2
itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
! nres_molec(molecule)=nres_molec(molecule)+1
+ read (card(19:19),'(a1)') sugar
+ isugar=sugarcode(sugar,ires)
+! if (ibeg.eq.1) then
+! istype(1)=isugar
+! else
+ istype(ires)=isugar
+! print *,"ires=",ires,istype(ires)
+! endif
+
endif
endif
else
.or. atom.eq."C4'" .or. atom.eq."O4'")) then
read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
!c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
- print *,ires,ishift,ishift1
+! print *,ires,ishift,ishift1
counter=counter+1
! iii=iii+1
! do j=1,3
counter=0
endif
! print *, "ATOM",atom(1:3)
- else if (atom(1:3).eq."C5'") then
+ else if (atom.eq."C5'") then
read (card(19:19),'(a1)') sugar
isugar=sugarcode(sugar,ires)
if (ibeg.eq.1) then
istype(1)=isugar
else
istype(ires)=isugar
+! print *,ires,istype(ires)
endif
if (unres_pdb) then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
if (itype(i,k).eq.ntyp1_molec(k)) then
if (itype(i+1,k).eq.ntyp1_molec(k)) then
if (itype(i+2,k).eq.0) then
- print *,"masz sieczke"
+! print *,"masz sieczke"
do j=1,5
if (itype(i+2,j).ne.0) then
itype(i+1,k)=0
e2(2)=1.0d0
e2(3)=0.0d0
endif !fail
- print *,i,'a tu?'
+! print *,i,'a tu?'
do j=1,3
c(j,i)=c(j,i-1)-1.9d0*e2(j)
enddo
allocate(c_temporary(3,2*nres))
allocate(itype_temporary(nres,5))
allocate(molnum(nres))
+ allocate(istype_temp(nres))
itype_temporary(:,:)=0
seqalingbegin=1
do k=1,5
enddo
itype_temporary(seqalingbegin,k)=itype(i,k)
+ istype_temp(seqalingbegin)=istype(i)
molnum(seqalingbegin)=k
seqalingbegin=seqalingbegin+1
endif
do k=1,5
do i=1,nres
itype(i,k)=itype_temporary(i,k)
+ istype(i)=istype_temp(i)
enddo
enddo
if (itype(1,1).eq.ntyp1) then
enddo
endif
- print *,seqalingbegin,nres
+! print *,seqalingbegin,nres
if(.not.allocated(vbld)) then
allocate(vbld(2*nres))
do i=1,2*nres
kkk=1
lll=0
cou=1
+ write (iout,*) "symetr", symetr
do i=1,nres
lll=lll+1
!c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
open (irotam_pdb,file=rotname_pdb,status='old',action='read')
#endif
#endif
+ call getenv_loc('SCPPAR_NUCL',scpname_nucl)
+#if defined(WINIFL) || defined(WINPGI)
+ open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
+#elif (defined G77)
+ open (iscpp_nucl,file=scpname_nucl,status='old')
+#else
+ open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
+#endif
+
#ifndef OLDSCP
!
! 8/9/01 In the newest version SCp interaction constants are read from a file
projects / unres4.git / commitdiff