From 04d0eb23476c67cc0cadc32bf09aa91c03e2dc15 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 9 Aug 2017 12:49:44 +0200 Subject: [PATCH] workink EVDW i EES for PP,SB,PSB- warning energies differ as corrections made for dummy atom --- PARAM/psb.parm | 5 + source/unres/control.F90 | 159 +++++++- source/unres/data/calc_data.f90 | 15 + source/unres/data/energy_data.f90 | 28 +- source/unres/data/io_units.f90 | 4 +- source/unres/data/names.f90 | 2 +- source/unres/energy.f90 | 763 ++++++++++++++++++++++++++++++++++++- source/unres/geometry.f90 | 64 +++- source/unres/io.f90 | 2 +- source/unres/io_config.f90 | 141 ++++++- 10 files changed, 1153 insertions(+), 30 deletions(-) create mode 100644 PARAM/psb.parm diff --git a/PARAM/psb.parm b/PARAM/psb.parm new file mode 100644 index 0000000..c9ec9f7 --- /dev/null +++ b/PARAM/psb.parm @@ -0,0 +1,5 @@ +-1.0 5.0 +-1.0 5.0 +-1.0 5.0 +-1.0 5.0 +-1.0 5.0 diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 86f23f9..b851d37 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -236,6 +236,7 @@ itordp_nucl= 130 ! ielep_nucl= 131 isidep_nucl=132 + iscpp_nucl=133 @@ -488,11 +489,14 @@ !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,& @@ -500,6 +504,11 @@ 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) @@ -662,6 +671,41 @@ 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 @@ -675,6 +719,12 @@ 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 @@ -707,6 +757,38 @@ 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) @@ -727,6 +809,31 @@ 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 @@ -751,6 +858,7 @@ 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) @@ -785,6 +893,47 @@ 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 @@ -933,10 +1082,10 @@ 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,& diff --git a/source/unres/data/calc_data.f90 b/source/unres/data/calc_data.f90 index ed6d106..7de07b0 100644 --- a/source/unres/data/calc_data.f90 +++ b/source/unres/data/calc_data.f90 @@ -12,3 +12,18 @@ 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 + diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index 844aa85..44e3d30 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -72,7 +72,7 @@ #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'/) @@ -93,20 +93,34 @@ 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) @@ -115,7 +129,19 @@ 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 diff --git a/source/unres/data/io_units.f90 b/source/unres/data/io_units.f90 index 8b79b35..5e9e9b7 100644 --- a/source/unres/data/io_units.f90 +++ b/source/unres/data/io_units.f90 @@ -16,7 +16,7 @@ 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 @@ -47,7 +47,7 @@ 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 diff --git a/source/unres/data/names.f90 b/source/unres/data/names.f90 index 97eb2d0..cb3f2f9 100644 --- a/source/unres/data/names.f90 +++ b/source/unres/data/names.f90 @@ -102,7 +102,7 @@ (/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 !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index fb4caac..051ba09 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -126,7 +126,8 @@ 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) @@ -549,6 +550,11 @@ 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 @@ -3188,10 +3194,10 @@ 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) @@ -19395,6 +19401,19 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -19567,6 +19586,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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)) @@ -19853,7 +19879,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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. & @@ -20013,7 +20039,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -20046,7 +20072,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) & @@ -20055,7 +20081,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -20100,7 +20126,726 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 6895472..e0b37d3 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -1617,6 +1617,68 @@ 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) @@ -3055,7 +3117,7 @@ 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 diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 9388f15..ee588d0 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -1187,7 +1187,7 @@ 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 diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 4c7f3a5..e6b7f0a 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -1683,7 +1683,7 @@ #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) @@ -1701,12 +1701,12 @@ !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) @@ -2055,6 +2055,7 @@ 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)) @@ -2164,6 +2165,7 @@ end select continue close (isidep) + !----------------------------------------------------------------------- ! Calculate the "working" parameters of SC interactions. @@ -2272,6 +2274,85 @@ 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 @@ -2344,6 +2425,20 @@ 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 ! @@ -2438,7 +2533,7 @@ ! 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 @@ -2475,6 +2570,7 @@ 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 @@ -2586,7 +2682,7 @@ 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 @@ -2606,6 +2702,15 @@ 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 @@ -2635,7 +2740,7 @@ .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 @@ -2653,13 +2758,14 @@ 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) @@ -2725,7 +2831,7 @@ 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 @@ -2749,7 +2855,7 @@ 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 @@ -2947,6 +3053,7 @@ 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 @@ -2958,6 +3065,7 @@ enddo itype_temporary(seqalingbegin,k)=itype(i,k) + istype_temp(seqalingbegin)=istype(i) molnum(seqalingbegin)=k seqalingbegin=seqalingbegin+1 endif @@ -2971,6 +3079,7 @@ 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 @@ -3008,7 +3117,7 @@ enddo endif - print *,seqalingbegin,nres +! print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -3087,6 +3196,7 @@ 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) @@ -4247,6 +4357,17 @@ 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 -- 1.7.9.5