! Get parameter filenames and open the parameter files.
call mygetenv('BONDPAR',bondname)
open (ibond,file=bondname,status='old')
+ call mygetenv('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old')
call mygetenv('THETPAR',thetname)
open (ithep,file=thetname,status='old')
call mygetenv('ROTPAR',rotname)
open (isidep,file=sidename,status='old')
call mygetenv('SIDEP',sidepname)
open (isidep1,file=sidepname,status="old")
+ call mygetenv('THETPAR_NUCL',thetname_nucl)
+ open (ithep_nucl,file=thetname_nucl,status='old')
+ call mygetenv('ROTPAR_NUCL',rotname_nucl)
+ open (irotam_nucl,file=rotname_nucl,status='old')
+ call mygetenv('TORPAR_NUCL',torname_nucl)
+ open (itorp_nucl,file=torname_nucl,status='old')
+ call mygetenv('TORDPAR_NUCL',tordname_nucl)
+ open (itordp_nucl,file=tordname_nucl,status='old')
+ call mygetenv('SIDEPAR_NUCL',sidename_nucl)
+ open (isidep_nucl,file=sidename_nucl,status='old')
+ call mygetenv('SIDEPAR_SCBASE',sidename_scbase)
+ open (isidep_scbase,file=sidename_scbase,status='old')
+ call mygetenv('PEPPAR_PEPBASE',pepname_pepbase)
+ open (isidep_pepbase,file=pepname_pepbase,status='old')
+ call mygetenv('SCPAR_PHOSPH',pepname_scpho)
+ open (isidep_scpho,file=pepname_scpho,status='old')
+ call mygetenv('PEPPAR_PHOSPH',pepname_peppho)
+ open (isidep_peppho,file=pepname_peppho,status='old')
+ call mygetenv('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
+ call mygetenv('TUBEPAR',tubename)
+ open (itube,file=tubename,status='old',action='read')
+ call mygetenv('IONPAR',ionname)
+ open (iion,file=ionname,status='old',action='read')
+
+ call mygetenv('SCPPAR_NUCL',scpname_nucl)
+ open (iscpp_nucl,file=scpname_nucl,status='old')
+ call mygetenv('IONPAR_NUCL',ionnuclname)
+ open (iionnucl,file=ionnuclname,status='old')
+ call mygetenv('IONPAR_TRAN',iontranname)
+ open (iiontran,file=iontranname,status='old',action='read')
+
+
#ifndef OLDSCP
!
! 8/9/01 In the newest version SCp interaction constants are read from a file
! Read molecular data.
!
use energy_data
- use geometry_data, only:nres,deg2rad,c,dc,nres_molec
- use control_data, only:iscode
+ use geometry_data, only:nres,deg2rad,c,dc,nres_molec,crefjlee,cref
+ use control_data, only:iscode,dyn_ss,pdbref,indpdb
use io_base, only:rescode
- use control, only:setup_var,init_int_table
+ use control, only:setup_var,init_int_table,hpb_partition
use geometry, only:alloc_geo_arrays
- use energy, only:alloc_ener_arrays
+ use energy, only:alloc_ener_arrays
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
endif
endif
-
nnt=1
nct=nres
if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
write(iout,*) 'NNT=',NNT,' NCT=',NCT
+ if (constr_homology.gt.0) then
+!c write (iout,*) "About to call read_constr_homology"
+!c call flush(iout)
+ call read_constr_homology
+!c write (iout,*) "Exit read_constr_homology"
+!c call flush(iout)
+ if (indpdb.gt.0 .or. pdbref) then
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=crefjlee(j,i)
+ cref(j,i,1)=crefjlee(j,i)
+ enddo
+ enddo
+ endif
+ endif
+#ifdef DEBUG
+ write (iout,*) "Array C"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
+ & (c(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Array Cref"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
+ & (cref(j,i+nres,1),j=1,3)
+ enddo
+#endif
call setup_var
call init_int_table
if (ns.gt.0) then
dhpb(i),ebr,forcon(i)
enddo
endif
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+ enddo
+ endif
write (iout,'(a)')
return
end subroutine molread
integer :: iparm,nkcctyp
!el real(kind=8) :: ip,mp
real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
- sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
+ sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, &
+ epspeptube,epssctube,sigmapeptube, &
+ sigmasctube,ssscale
+
real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
res1
integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
character*3 string
-
+
!
! Body
!
vbl=3.8D0
vblinv=1.0D0/vbl
vblinv2=vblinv*vblinv
+ itime_mat=0
#ifndef CLUSTER
call card_concat(controlcard,.true.)
wname(4)="WCORRH"
!el
+ call reada(controlcard,"D0CM",d0cm,3.78d0)
+ call reada(controlcard,"AKCM",akcm,15.1d0)
+ call reada(controlcard,"AKTH",akth,11.0d0)
+ call reada(controlcard,"AKCT",akct,12.0d0)
+ call reada(controlcard,"V1SS",v1ss,-1.08d0)
+ call reada(controlcard,"V2SS",v2ss,7.61d0)
+ call reada(controlcard,"V3SS",v3ss,13.7d0)
+ call reada(controlcard,"EBR",ebr,-5.50D0)
+ call reada(controlcard,"ATRISS",atriss,0.301D0)
+ call reada(controlcard,"BTRISS",btriss,0.021D0)
+ call reada(controlcard,"CTRISS",ctriss,1.001D0)
+ call reada(controlcard,"DTRISS",dtriss,1.001D0)
+ call reada(controlcard,"SSSCALE",ssscale,1.0D0)
+ dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+ call reada(controlcard,"LIPSCALE",lipscale,1.0D0)
allocate(ww(max_eneW))
do i=1,n_eneW
key = wname(i)(:ilen(wname(i)))
wscloc=ww(12)
wtor=ww(13)
wtor_d=ww(14)
+ wstrain=ww(15)
wvdwpp=ww(16)
wbond=ww(18)
wsccor=ww(19)
wcatcat=ww(42)
wcatprot=ww(41)
+ wcorr3_nucl=ww(38)
+ wcorr_nucl=ww(37)
+ wtor_d_nucl=ww(36)
+ wtor_nucl=ww(35)
+ wsbloc=ww(34)
+ wang_nucl=ww(33)
+ wbond_nucl=ww(32)
+ welsb=ww(31)
+ wvdwsb=ww(30)
+ welpsb=ww(29)
+ wvdwpsb=ww(28)
+ welpp=ww(27)
+ wvdwpp_nucl=ww(26)
+ wscbase=ww(46)
+ wpepbase=ww(47)
+ wscpho=ww(48)
+ wpeppho=ww(49)
+ wcatnucl=ww(50)
+ wcat_tran=ww(56)
+! print *,"KURWA",ww(48)
+! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "
+! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
+! "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",&
+! "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",&
+! "WCATPROT ","WCATCAT
+! +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+! +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+! +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+! +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
endif
!
weights(12)=wscloc
weights(13)=wtor
weights(14)=wtor_d
- weights(15)=0 !wstrain !
- weights(16)=0 !wvdwpp !
+ weights(15)=wstrain !0
+ weights(16)=wvdwpp !
weights(17)=wbond
weights(18)=0 !scal14 !
weights(21)=wsccor
weights(42)=wcatprot
weights(41)=wcatcat
+ weights(26)= wvdwpp_nucl
+
+ weights(27) =welpp
+ weights(28) =wvdwpsb
+ weights(29) =welpsb
+ weights(30) =wvdwsb
+ weights(31) =welsb
+ weights(32) =wbond_nucl
+ weights(33) =wang_nucl
+ weights(34) =wsbloc
+ weights(35) =wtor_nucl
+ weights(36) =wtor_d_nucl
+ weights(37) =wcorr_nucl
+ weights(38) =wcorr3_nucl
+ weights(41) =wcatcat
+ weights(42) =wcatprot
+ weights(46) =wscbase
+ weights(47)= wpepbase
+ weights(48) =wscpho
+ weights(49) =wpeppho
+ weights(50) =wcatnucl
+ weights(56)=wcat_tran
! el--------
call card_concat(controlcard,.false.)
write (iout,*) "Parameter set:",iparm
write (iout,*) "Energy-term weights:"
do i=1,n_eneW
- write (iout,'(a16,f10.5)') wname(i),ww(i)
+ write (iout,'(i3,a16,f10.5)') i,wname(i),ww(i)
enddo
write (iout,*) "Sidechain potential file : ",&
sidename_t(:ilen(sidename_t))
allocate(nbondterm(ntyp)) !(ntyp)
allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+ allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
+ allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+ allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+
!el allocate(msc(ntyp+1)) !(ntyp+1)
!el allocate(isc(ntyp+1)) !(ntyp+1)
!el allocate(restok(ntyp+1)) !(ntyp+1)
enddo
enddo
endif
- if (.not. allocated(msc)) allocate(msc(ntyp1,5))
- if (.not. allocated(restok)) allocate(restok(ntyp1,5))
+ if (.not. allocated(msc)) allocate(msc(-ntyp1:ntyp1,5))
+ if (.not. allocated(restok)) allocate(restok(-ntyp1:ntyp1,5))
+ if (oldion.eq.1) then
do i=1,ntyp_molec(5)
read(iion,*) msc(i,5),restok(i,5)
read (iion,*) (catprm(i,k),i=1,ncatprotparm)
enddo
print *, catprm
+ endif
+ allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
+ do i=1,ntyp_molec(5)
+ do j=1,ntyp_molec(2)
+ write(iout,*) i,j
+ read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
+ enddo
+ enddo
+ write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
+ "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
+ "b3 ","b4 ","chis1 "
+ do i=1,ntyp_molec(5)
+ do j=1,ntyp_molec(2)
+ write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
+ restyp(i,5),"-",restyp(j,2)
+ enddo
+ enddo
+
+ read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
+ do i=1,ntyp_molec(2)
+ nbondterm_nucl(i)=1
+ read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
+! dsc(i) = vbldsc0_nucl(1,i)
+! if (i.eq.10) then
+! dsc_inv(i)=0.0D0
+! else
+! dsc_inv(i)=1.0D0/dsc(i)
+! endif
+ enddo
!----------------------------------------------------
ENDIF ! TOR_MODE
#endif
+!--------------- Reading theta parameters for nucleic acid-------
+ read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
+ ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
+ nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
+ allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ allocate(aa0thet_nucl(nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxtheterm,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+
+ read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
+
+ aa0thet_nucl(:,:,:)=0.0d0
+ aathet_nucl(:,:,:,:)=0.0d0
+ bbthet_nucl(:,:,:,:,:)=0.0d0
+ ccthet_nucl(:,:,:,:,:)=0.0d0
+ ddthet_nucl(:,:,:,:,:)=0.0d0
+ eethet_nucl(:,:,:,:,:)=0.0d0
+ ffthet_nucl(:,:,:,:,:,:)=0.0d0
+ ggthet_nucl(:,:,:,:,:,:)=0.0d0
+
+ do i=1,nthetyp_nucl
+ do j=1,nthetyp_nucl
+ do k=1,nthetyp_nucl
+ read (ithep_nucl,'(3a)') t1,t2,t3
+ read (ithep_nucl,*) aa0thet_nucl(i,j,k)
+ read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
+ read (ithep_nucl,*) &
+ (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
+ read (ithep_nucl,*) &
+ (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
+ ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
+ llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
+ enddo
+ enddo
+ enddo
+
+
!-------------------------------------------
allocate(nlob(ntyp1)) !(ntyp1)
allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
enddo
#endif
close(irotam)
+!---------reading nucleic acid parameters for rotamers-------------------
+ allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
+ do i=1,ntyp_molec(2)
+ read (irotam_nucl,*)
+ do j=1,9
+ read(irotam_nucl,*) sc_parmin_nucl(j,i)
+ enddo
+ enddo
+ close(irotam_nucl)
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) "Base rotamer parameters"
+ do i=1,ntyp_molec(2)
+ write (iout,'(a)') restyp(i,2)
+ write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
+ enddo
+ endif
+
read (ifourier,*) nloctyp
!el write(iout,*)"nloctyp",nloctyp
enddo
enddo
endif
+#ifndef NEWCORR
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ enddo
+#endif
ELSE IF (TOR_MODE.eq.1) THEN
!C read valence-torsional parameters
allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
+ allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
+ dcavtub(ntyp1))
+ allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
+ tubetranene(ntyp1))
do i=1,ntyp1
do j=1,ntyp1
aa_aq(i,j)=0.0D0
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
+ sigmapeptube=sigmapeptube**6
+ sigeps=dsign(1.0D0,epspeptube)
+ epspeptube=dabs(epspeptube)
+ pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
+ pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
+ write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
+ do i=1,ntyp
+ read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
+ ccavtub(i),dcavtub(i),tubetranene(i)
+ sigmasctube=sigmasctube**6
+ sigeps=dsign(1.0D0,epssctube)
+ epssctube=dabs(epssctube)
+ sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
+ sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
+ write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
+ enddo
+!-----------------READING SC BASE POTENTIALS-----------------------------
+ allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
+ allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
+ allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
+ allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
+ allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
+ allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
+ allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
+ allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
+ allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
+ allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
+ allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+
+ do i=1,ntyp_molec(1)
+ do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+ write (*,*) "Im in ", i, " ", j
+ read(isidep_scbase,*) &
+ eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+ chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+ write(*,*) "eps",eps_scbase(i,j)
+ read(isidep_scbase,*) &
+ (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+ chis_scbase(i,j,1),chis_scbase(i,j,2)
+ read(isidep_scbase,*) &
+ dhead_scbasei(i,j), &
+ dhead_scbasej(i,j), &
+ rborn_scbasei(i,j),rborn_scbasej(i,j)
+ read(isidep_scbase,*) &
+ (wdipdip_scbase(k,i,j),k=1,3), &
+ (wqdip_scbase(k,i,j),k=1,2)
+ read(isidep_scbase,*) &
+ alphapol_scbase(i,j), &
+ epsintab_scbase(i,j)
+ END DO
+ END DO
+ allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+ do i=1,ntyp_molec(1)
+ do j=1,ntyp_molec(2)-1
+ epsij=eps_scbase(i,j)
+ rrij=sigma_scbase(i,j)
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_scbase(i,j)=epsij*rrij*rrij
+ bb_scbase(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+!-----------------READING PEP BASE POTENTIALS-------------------
+ allocate(eps_pepbase(ntyp_molec(2)))
+ allocate(sigma_pepbase(ntyp_molec(2)))
+ allocate(chi_pepbase(ntyp_molec(2),2))
+ allocate(chipp_pepbase(ntyp_molec(2),2))
+ allocate(alphasur_pepbase(4,ntyp_molec(2)))
+ allocate(sigmap1_pepbase(ntyp_molec(2)))
+ allocate(sigmap2_pepbase(ntyp_molec(2)))
+ allocate(chis_pepbase(ntyp_molec(2),2))
+ allocate(wdipdip_pepbase(3,ntyp_molec(2)))
+
+
+ do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+ write (*,*) "Im in ", i, " ", j
+ read(isidep_pepbase,*) &
+ eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+ chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+ write(*,*) "eps",eps_pepbase(j)
+ read(isidep_pepbase,*) &
+ (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+ chis_pepbase(j,1),chis_pepbase(j,2)
+ read(isidep_pepbase,*) &
+ (wdipdip_pepbase(k,j),k=1,3)
+ END DO
+ allocate(aa_pepbase(ntyp_molec(2)))
+ allocate(bb_pepbase(ntyp_molec(2)))
+
+ do j=1,ntyp_molec(2)-1
+ epsij=eps_pepbase(j)
+ rrij=sigma_pepbase(j)
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_pepbase(j)=epsij*rrij*rrij
+ bb_pepbase(j)=-sigeps*epsij*rrij
+ enddo
+!--------------READING SC PHOSPHATE-------------------------------------
+!--------------READING SC PHOSPHATE-------------------------------------
+ allocate(eps_scpho(ntyp_molec(1)))
+ allocate(sigma_scpho(ntyp_molec(1)))
+ allocate(chi_scpho(ntyp_molec(1),2))
+ allocate(chipp_scpho(ntyp_molec(1),2))
+ allocate(alphasur_scpho(4,ntyp_molec(1)))
+ allocate(sigmap1_scpho(ntyp_molec(1)))
+ allocate(sigmap2_scpho(ntyp_molec(1)))
+ allocate(chis_scpho(ntyp_molec(1),2))
+ allocate(wqq_scpho(ntyp_molec(1)))
+ allocate(wqdip_scpho(2,ntyp_molec(1)))
+ allocate(alphapol_scpho(ntyp_molec(1)))
+ allocate(epsintab_scpho(ntyp_molec(1)))
+ allocate(dhead_scphoi(ntyp_molec(1)))
+ allocate(rborn_scphoi(ntyp_molec(1)))
+ allocate(rborn_scphoj(ntyp_molec(1)))
+ allocate(alphi_scpho(ntyp_molec(1)))
+
+
+! j=1
+ do j=1,ntyp_molec(1) ! without U then we will take T for U
+ write (*,*) "Im in scpho ", i, " ", j
+ read(isidep_scpho,*) &
+ eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
+ chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
+ write(*,*) "eps",eps_scpho(j)
+ read(isidep_scpho,*) &
+ (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
+ chis_scpho(j,1),chis_scpho(j,2)
+ read(isidep_scpho,*) &
+ (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
+ read(isidep_scpho,*) &
+ epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
+ alphi_scpho(j)
+
+ END DO
+ allocate(aa_scpho(ntyp_molec(1)))
+ allocate(bb_scpho(ntyp_molec(1)))
+ do j=1,ntyp_molec(1)
+ epsij=eps_scpho(j)
+ rrij=sigma_scpho(j)
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_scpho(j)=epsij*rrij*rrij
+ bb_scpho(j)=-sigeps*epsij*rrij
+ enddo
+
+
+ read(isidep_peppho,*) &
+ eps_peppho,sigma_peppho
+ read(isidep_peppho,*) &
+ (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
+ read(isidep_peppho,*) &
+ (wqdip_peppho(k),k=1,2)
+
+ epsij=eps_peppho
+ rrij=sigma_peppho
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_peppho=epsij*rrij*rrij
+ bb_peppho=-sigeps*epsij*rrij
+
allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
do i=1,ntyp
bad(i,j)=0.0D0
enddo
enddo
+! Ions by Aga
+
+ allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+ allocate(alphapolcat2(ntyp,ntyp))
+ allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+ allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+ allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
+ allocate(epsintabcat(ntyp,ntyp))
+ allocate(dtailcat(2,ntyp,ntyp))
+ allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
+ allocate(wqdipcat(2,ntyp,ntyp))
+ allocate(wstatecat(4,ntyp,ntyp))
+ allocate(dheadcat(2,2,ntyp,ntyp))
+ allocate(nstatecat(ntyp,ntyp))
+ allocate(debaykapcat(ntyp,ntyp))
+
+ if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
+ if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+
+ allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
+! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
+ if (oldion.eq.0) then
+ if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
+ allocate(icharge(1:ntyp1))
+ read(iion,*) (icharge(i),i=1,ntyp)
+ else
+ read(iion,*) ijunk
+ endif
+
+ do i=-ntyp_molec(5),ntyp_molec(5)
+ read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
+ print *,msc(i,5),restok(i,5)
+ enddo
+ ip(5)=0.2
+!DIR$ NOUNROLL
+ do j=1,ntyp_molec(5)-1
+ do i=1,ntyp
+! do j=1,ntyp_molec(5)
+! write (*,*) "Im in ALAB", i, " ", j
+ read(iion,*) &
+ epscat(i,j),sigmacat(i,j), &
+! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
+ chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+
+ (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+! chiscat(i,j),chiscat(j,i), &
+ chis1cat(i,j),chis2cat(i,j), &
+
+ nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
+ dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+ dtailcat(1,i,j),dtailcat(2,i,j), &
+ epsheadcat(i,j),sig0headcat(i,j), &
+!wdipcat = w1 , w2
+! rborncat(i,j),rborncat(j,i),&
+ rborn1cat(i,j),rborn2cat(i,j),&
+ (wqdipcat(k,i,j),k=1,2), &
+ alphapolcat(i,j),alphapolcat2(j,i), &
+ (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
+! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
+ END DO
+ END DO
+ allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)-1 !without zinc
+ epsij=epscat(i,j)
+ rrij=sigmacat(i,j)
+ rrij=rrij**expon
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_aq_cat(i,j)=epsij*rrij*rrij
+ bb_aq_cat(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)
+ if (i.eq.10) then
+ write (iout,*) 'i= ', i, ' j= ', j
+ write (iout,*) 'epsi0= ', epscat(i,j)
+ write (iout,*) 'sigma0= ', sigmacat(i,j)
+ write (iout,*) 'chi1= ', chi1cat(i,j)
+ write (iout,*) 'chi1= ', chi2cat(i,j)
+ write (iout,*) 'chip1= ', chipp1cat(1,j)
+ write (iout,*) 'chip2= ', chipp2cat(1,j)
+ write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
+ write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
+ write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
+ write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
+ write (iout,*) 'sig1= ', sigmap1cat(1,j)
+ write (iout,*) 'sig2= ', sigmap2cat(1,j)
+ write (iout,*) 'chis1= ', chis1cat(1,j)
+ write (iout,*) 'chis1= ', chis2cat(1,j)
+ write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
+ write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
+ write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
+ write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
+ write (iout,*) 'a1= ', rborn1cat(i,j)
+ write (iout,*) 'a2= ', rborn2cat(i,j)
+ write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
+ write (iout,*) 'alphapol1= ', alphapolcat(1,j)
+ write (iout,*) 'alphapol2= ', alphapolcat(j,1)
+ write (iout,*) 'w1= ', wqdipcat(1,i,j)
+ write (iout,*) 'w2= ', wqdipcat(2,i,j)
+ write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
+ endif
+
+ If ((i.eq.1).and.(j.eq.27)) then
+ write (iout,*) 'SEP'
+ Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+ Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+ endif
+
+ enddo
+ enddo
+
+ endif
+! here two denotes the Zn2+ and Cu2+
+ write(iout,*) "before TRANPARM"
+ allocate(aomicattr(0:3,2))
+ allocate(athetacattran(0:6,5,2))
+ allocate(agamacattran(3,5,2))
+ allocate(acatshiftdsc(5,2))
+ allocate(bcatshiftdsc(5,2))
+ allocate(demorsecat(5,2))
+ allocate(alphamorsecat(5,2))
+ allocate(x0catleft(5,2))
+ allocate(x0catright(5,2))
+ allocate(x0cattrans(5,2))
+ allocate(ntrantyp(2))
+ do i=1,1 ! currently only Zn2+
+
+ read(iiontran,*) ntrantyp(i)
+!ntrantyp=4
+!| ao0 ao1 ao2 ao3
+!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
+!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
+!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+ read(iiontran,*) (aomicattr(j,i),j=0,3)
+ do j=1,ntrantyp(i)
+ read (iiontran,*) (agamacattran(k,j,i),k=1,3),&
+ (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+ demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+ x0cattrans(j,i)
+ enddo
+ enddo
#ifdef CLUSTER
!
! Define the SC-p interaction constants
enddo
enddo
#endif
+ allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ read (itorp_nucl,*) ntortyp_nucl
+! 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)
+
+ allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+ allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
+ allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
+ allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
+ allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el---------------------------
+ nterm_nucl(:,:)=0
+ nlor_nucl(:,:)=0
+!el--------------------
+ read (itorp_nucl,*) &
+ (itortyp_nucl(i),i=1,ntyp_molec(2))
+! print *,itortyp_nucl(:)
+!c write (iout,*) 'ntortyp',ntortyp
+ do i=1,ntortyp_nucl
+ do j=1,ntortyp_nucl
+ read (itorp_nucl,*) 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)
+ read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
+ v0ij=v0ij+si*v1_nucl(k,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_nucl(i,j)
+ read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
+ vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
+ v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
+ enddo
+ v0_nucl(i,j)=v0ij
+ enddo
+ enddo
+
!elwrite(iout,*) "parmread kontrol before oldscp"
!
enddo
endif
#endif
+ allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
+
+ do i=1,ntyp_molec(2)
+ read (iscpp_nucl,*) 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
!
- ebr=-5.50D0
+! ebr=-5.50D0
!
! Old arbitrary potential - commented out.
!
! energy surface of diethyl disulfide.
! A. Liwo and U. Kozlowska, 11/24/03
!
- D0CM = 3.78d0
- AKCM = 15.1d0
- AKTH = 11.0d0
- AKCT = 12.0d0
- V1SS =-1.08d0
- V2SS = 7.61d0
- V3SS = 13.7d0
+! D0CM = 3.78d0
+! AKCM = 15.1d0
+! AKTH = 11.0d0
+! AKCT = 12.0d0
+! V1SS =-1.08d0
+! V2SS = 7.61d0
+! V3SS = 13.7d0
+#ifndef CLUSTER
+ if (dyn_ss) then
+ ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
+ Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
+ akcm=akcm/wsc*ssscale
+ akth=akth/wsc*ssscale
+ akct=akct/wsc*ssscale
+ v1ss=v1ss/wsc*ssscale
+ v2ss=v2ss/wsc*ssscale
+ v3ss=v3ss/wsc*ssscale
+ else
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+#endif
if (lprint) then
write (iout,'(/a)') "Disulfide bridge parameters:"
write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
subroutine read_general_data(*)
use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
- scelemode,TUBEmode,tor_mode
+ scelemode,TUBEmode,tor_mode,energy_dec,r_cut_ang,r_cut_mart,&
+ rlamb_mart
- use energy_data, only:distchainmax,tubeR0,tubecenter
+ use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss,constr_homology
use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
- bordtubebot,tubebufthick,buftubebot,buftubetop
+ bordtubebot,tubebufthick,buftubebot,buftubetop,buflipbot, bufliptop,bordlipbot,bordliptop, &
+ lipbufthick,lipthick
! implicit none
! include "DIMENSIONS"
! include "DIMENSIONS.ZSCOPT"
! include "COMMON.FREE"
! include "COMMON.CONTROL"
! include "COMMON.ENERGIES"
+
character(len=800) :: controlcard
integer :: i,j,k,ii,n_ene_found
integer :: ind,itype1,itype2,itypf,itypsc,itypp
call reada(controlcard,'BOXX',boxxsize,100.0d0)
call reada(controlcard,'BOXY',boxysize,100.0d0)
call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+ write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick) &
+ write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif !lipthick.gt.0
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+
+ energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
call readi(controlcard,"SCELEMODE",scelemode,0)
+ call readi(controlcard,"OLDION",oldion,0)
+ dyn_ss=(index(controlcard,'DYN_SS').gt.0)
print *,"SCELE",scelemode
call readi(controlcard,'TORMODE',tor_mode,0)
!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write(iout,*) "torsional and valence angle mode",tor_mode
call readi(controlcard,'TUBEMOD',tubemode,0)
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+!c if (constr_homology) tole=dmax1(tole,1.5d0)
+ write (iout,*) "with_homology_constr ",with_dihed_constr, &
+ " CONSTR_HOMOLOGY",constr_homology
+! read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
+! out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
+! out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
if (TUBEmode.gt.0) then
call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
buftubetop=bordtubetop-tubebufthick
endif
ions=index(controlcard,"IONS").gt.0
- call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+ call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
write(iout,*) "R_CUT_ELE=",r_cut_ele
+ call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
+ call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
+ call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
call readi(controlcard,'SYM',symetr,1)
do k=1,3
cref(k,j+i,kkk)=cref(k,j,kkk)
enddo
+ write(iout,*) "J",j,"J+I",j+i
phi_ref(j+i)=phi_ref(j)
theta_ref(j+i)=theta_ref(j)
alph_ref(j+i)=alph_ref(j)
!--------------------------------------------------------------------------------
subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
- use geometry_data, only:nres,c
+ use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize
use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
+ use energy, only:boxshift
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
! include 'COMMON.HEADER'
! include 'COMMON.SBRIDGE'
character(len=50) :: tytul
- character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
- 'D','E','F','G','H','I','J'/),shape(chainid))
+ character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
+ 'D','E','F','G','H','I','J','K','L','M','N','O',&
+ 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
integer,dimension(nres) :: ica !(maxres)
- real(kind=8) :: temp,efree,etot,entropy,rmsdev
+ real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj
integer :: ii,i,j,iti,ires,iatom,ichain,mnum
write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
ii,temp,rmsdev
do i=nnt,nct
mnum=molnum(i)
iti=itype(i,mnum)
- if (iti.eq.ntyp1) then
+ if (iti.eq.ntyp1_molec(mnum)) then
+ if (itype(i-1,molnum(i-1)).eq.ntyp1_molec(mnum)) then
ichain=ichain+1
ires=0
write (ipdb,'(a)') 'TER'
+ endif
else
ires=ires+1
iatom=iatom+1
ica(i)=iatom
+ if (mnum.ne.5) then
write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,i),j=1,3)
- if (iti.ne.10) then
+ else
+ xj=boxshift(c(1,i)-c(1,2),boxxsize)
+ yj=boxshift(c(2,i)-c(2,2),boxysize)
+ zj=boxshift(c(3,i)-c(3,2),boxzsize)
+ write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
+ ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj
+ endif
+ if ((iti.ne.10).and.(mnum.ne.5)) then
iatom=iatom+1
write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,nres+i),j=1,3)
write (ipdb,'(a)') 'TER'
do i=nnt,nct-1
mnum=molnum(i)
- if (itype(i,mnum).eq.ntyp1) cycle
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
write (ipdb,30) ica(i),ica(i+1)
else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
return
end subroutine pdboutW
#endif
+
+ subroutine read_constr_homology
+ use energy_data
+ use control, only:init_int_table,homology_partition
+ use MD_data, only:iset
+ use geometry_data !only:nres,deg2rad,c,dc,nres_molec,crefjlee,cref
+ use MPI_data, only:kolor
+ use io_config, only:readpdb_template,readpdb
+! implicit none
+! include 'DIMENSIONS'
+!#ifdef MPI
+! include 'mpif.h'
+!#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.CONTROL'
+! include 'COMMON.HOMOLOGY'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+! include 'COMMON.QRESTR'
+! include 'COMMON.GEO'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.VAR'
+!
+
+! double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+! & dist_cut
+! common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+! & sigma_odl_temp(maxres,maxres,max_template)
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ integer :: ki,i,ii,j,k,l
+ integer, dimension (:), allocatable :: ii_in_use
+ integer :: i_tmp,idomain_tmp,&
+ irec,ik,iistart,nres_temp
+! integer :: iset
+! external :: ilen
+ logical :: liiflag,lfirst
+ integer :: i01,i10
+!
+! FP - Nov. 2014 Temporary specifications for new vars
+!
+ real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
+ rescore3_tmp, dist_cut
+ real(kind=8), dimension (:,:),allocatable :: rescore
+ real(kind=8), dimension (:,:),allocatable :: rescore2
+ real(kind=8), dimension (:,:),allocatable :: rescore3
+ real(kind=8) :: distal
+ character*24 tpl_k_rescore
+ character*256 pdbfile
+
+! -----------------------------------------------------------------
+! Reading multiple PDB ref structures and calculation of retraints
+! not using pre-computed ones stored in files model_ki_{dist,angle}
+! FP (Nov., 2014)
+! -----------------------------------------------------------------
+!
+!
+! Alternative: reading from input
+ call card_concat(controlcard,.true.)
+ call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
+ call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
+ call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
+ call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
+ call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
+ call readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
+ start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
+! if(.not.read2sigma.and.start_from_model) then
+! if(me1.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
+! write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
+! start_from_model=.false.
+! iranconf=(indpdb.le.0)
+! endif
+! if(start_from_model)&! .and. (me1.eq.king .or. .not. out1file))&
+! write(iout,*) 'START_FROM_MODELS is ON'
+! if(start_from_model .and. rest) then
+! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+! write(iout,*) 'START_FROM_MODELS is OFF'
+! write(iout,*) 'remove restart keyword from input'
+! endif
+! endif
+! if (start_from_model) nmodel_start=constr_homology
+ if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
+ if (homol_nset.gt.1)then
+ call card_concat(controlcard,.true.)
+ read(controlcard,*) (waga_homology(i),i=1,homol_nset)
+! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+! write(iout,*) "iset homology_weight "
+! do i=1,homol_nset
+ ! write(iout,*) i,waga_homology(i)
+ ! enddo
+ ! endif
+ iset=mod(kolor,homol_nset)+1
+ else
+ iset=1
+ waga_homology(1)=1.0
+ endif
+
+!d write (iout,*) "nnt",nnt," nct",nct
+!d call flush(iout)
+ if (.false.) then
+ print *,"here klapaciuj"
+! if (read_homol_frag) then
+! call read_klapaucjusz
+ else
+
+ lim_odl=0
+ lim_dih=0
+!
+! write(iout,*) 'nnt=',nnt,'nct=',nct
+!
+! do i = nnt,nct
+! do k=1,constr_homology
+! idomain(k,i)=0
+! enddo
+! enddo
+! idomain=0
+
+! ii=0
+! do i = nnt,nct-2
+! do j=i+2,nct
+! ii=ii+1
+! ii_in_use(ii)=0
+! enddo
+! enddo
+ ii_in_use=0
+ if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
+ if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
+
+ do k=1,constr_homology
+
+ read(inp,'(a)') pdbfile
+ pdbfiles_chomo(k)=pdbfile
+! if(me.eq.king .or. .not. out1file) &
+ write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
+ pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
+ pdbfile(:ilen(pdbfile))
+ stop
+ 34 continue
+! print *,'Begin reading pdb data'
+!
+! Files containing res sim or local scores (former containing sigmas)
+!
+
+ write(kic2,'(bz,i2.2)') k
+
+ tpl_k_rescore="template"//kic2//".sco"
+ write(iout,*) "tpl_k_rescore=",tpl_k_rescore
+! unres_pdb=.false.
+ nres_temp=nres
+ write(iout,*) "read2sigma",read2sigma
+
+ if (read2sigma) then
+ call readpdb_template(k)
+ else
+ call readpdb
+ endif
+ write(iout,*) "after readpdb"
+ if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
+ nres_chomo(k)=nres
+ nres=nres_temp
+ if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
+ if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
+ if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
+ if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
+ if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
+ if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
+ if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
+ if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
+ if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
+ if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
+ if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
+ if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
+ if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
+ if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
+! if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
+ if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
+ if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
+ if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
+ if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
+! if(.not.allocated(distance)) allocate(distance(constr_homology))
+! if(.not.allocated(distancek)) allocate(distancek(constr_homology))
+
+
+!
+! Distance restraints
+!
+! ... --> odl(k,ii)
+! Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres_chomo(k)
+ do j=1,3
+ c(j,i)=cref(j,i,1)
+! write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+!
+! From read_dist_constr (commented out 25/11/2014 <-> res sim)
+!
+! write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
+ open (ientin,file=tpl_k_rescore,status='old')
+ if (nnt.gt.1) rescore(k,1)=0.0d0
+ do irec=nnt,nct ! loop for reading res sim
+ if (read2sigma) then
+ read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
+ rescore3_tmp,idomain_tmp
+ i_tmp=i_tmp+nnt-1
+ write (*,*) "i_tmp", i_tmp,nnt
+ idomain(k,i_tmp)=idomain_tmp
+ rescore(k,i_tmp)=rescore_tmp
+ rescore2(k,i_tmp)=rescore2_tmp
+ rescore3(k,i_tmp)=rescore3_tmp
+! if (.not. out1file .or. me.eq.king)&
+! write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
+! i_tmp,rescore2_tmp,rescore_tmp,&
+! rescore3_tmp,idomain_tmp
+ else
+ idomain(k,irec)=1
+ read (ientin,*,end=1401) rescore_tmp
+
+! rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
+ rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
+! write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+ endif
+ enddo
+ 1401 continue
+ close (ientin)
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+
+ x12=c(1,i)-c(1,j)
+ y12=c(2,i)-c(2,j)
+ z12=c(3,i)-c(3,j)
+ distal=dsqrt(x12*x12+y12*y12+z12*z12)
+! write (iout,*) k,i,j,distal,dist2_cut
+
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
+ .and. distal.le.dist2_cut ) then
+
+ ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
+! write (iout,*) "k",k
+! write (iout,*) "i",i," j",j," constr_homology",
+! & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ odl(k,ii)=distal
+ if (read2sigma) then
+ sigma_odl(k,ii)=0
+ do ik=i,j
+ sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+ enddo
+ sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+ if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
+ sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+ else
+ if (odl(k,ii).le.dist_cut) then
+ sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
+ else
+#ifdef OLDSIGMA
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
+ dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
+ dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+#endif
+ endif
+ endif
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
+ else
+! ii=ii+1
+! l_homo(k,ii)=.false.
+ endif
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+! write (iout,*) "Distance restraints set"
+! call flush(iout)
+!
+! Theta, dihedral and SC retraints
+!
+ if (waga_angle.gt.0.0d0) then
+! open (ientin,file=tpl_k_sigma_dih,status='old')
+! do irec=1,maxres-3 ! loop for reading sigma_dih
+! read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+! if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+! sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+! & sigma_dih(k,i+nnt-1)
+! enddo
+!1402 continue
+! close (ientin)
+ do i = nnt+3,nct
+ if (idomain(k,i).eq.0) then
+ sigma_dih(k,i)=0.0
+ cycle
+ endif
+ dih(k,i)=phiref(i) ! right?
+! read (ientin,*) sigma_dih(k,i) ! original variant
+! write (iout,*) "dih(",k,i,") =",dih(k,i)
+! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+! & "rescore(",k,i-1,") =",rescore(k,i-1),
+! & "rescore(",k,i-2,") =",rescore(k,i-2),
+! & "rescore(",k,i-3,") =",rescore(k,i-3)
+
+ sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
+ rescore(k,i-2)+rescore(k,i-3))/4.0
+! if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
+! write (iout,*) "Raw sigmas for dihedral angle restraints"
+! write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
+! sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+! rescore(k,i-2)*rescore(k,i-3) ! right expression ?
+! Instead of res sim other local measure of b/b str reliability possible
+ if (sigma_dih(k,i).ne.0) &
+ sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+! sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
+ enddo
+ lim_dih=nct-nnt-2
+ endif
+! write (iout,*) "Dihedral angle restraints set"
+! call flush(iout)
+
+ if (waga_theta.gt.0.0d0) then
+! open (ientin,file=tpl_k_sigma_theta,status='old')
+! do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+! read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+! sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+! & sigma_theta(k,i+nnt-1)
+! enddo
+!1403 continue
+! close (ientin)
+
+ do i = nnt+2,nct ! right? without parallel.
+! do i = i=1,nres ! alternative for bounds acc to readpdb?
+! do i=ithet_start,ithet_end ! with FG parallel.
+ if (idomain(k,i).eq.0) then
+ sigma_theta(k,i)=0.0
+ cycle
+ endif
+ thetatpl(k,i)=thetaref(i)
+! write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+! & "rescore(",k,i-1,") =",rescore(k,i-1),
+! & "rescore(",k,i-2,") =",rescore(k,i-2)
+! read (ientin,*) sigma_theta(k,i) ! 1st variant
+ sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
+ rescore(k,i-2))/3.0
+! if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
+ if (sigma_theta(k,i).ne.0) &
+ sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+
+! sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+! rescore(k,i-2) ! right expression ?
+! sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
+ enddo
+ endif
+! write (iout,*) "Angle restraints set"
+! call flush(iout)
+
+ if (waga_d.gt.0.0d0) then
+! open (ientin,file=tpl_k_sigma_d,status='old')
+! do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+! read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+! sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+! & sigma_d(k,i+nnt-1)
+! enddo
+!1404 continue
+
+ do i = nnt,nct ! right? without parallel.
+! do i=2,nres-1 ! alternative for bounds acc to readpdb?
+! do i=loc_start,loc_end ! with FG parallel.
+ if (itype(i,1).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+ sigma_d(k,i)=0.0
+ cycle
+ endif
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+! write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
+! write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
+! write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
+! write(iout,*) "rescore(",k,i,") =",rescore(k,i)
+ sigma_d(k,i)=rescore3(k,i) ! right expression ?
+ if (sigma_d(k,i).ne.0) &
+ sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+
+! sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
+! sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
+! read (ientin,*) sigma_d(k,i) ! 1st variant
+ enddo
+ endif
+ enddo
+! write (iout,*) "SC restraints set"
+! call flush(iout)
+!
+! remove distance restraints not used in any model from the list
+! shift data in all arrays
+!
+! write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ liiflag=.true.
+ lfirst=.true.
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+! & .and. distal.le.dist2_cut ) then
+! write (iout,*) "i",i," j",j," ii",ii
+! call flush(iout)
+ if (ii_in_use(ii).eq.0.and.liiflag.or. &
+ ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
+ liiflag=.false.
+ i10=ii
+ if (lfirst) then
+ lfirst=.false.
+ iistart=ii
+ else
+ if(i10.eq.lim_odl) i10=i10+1
+ do ki=0,i10-i01-1
+ ires_homo(iistart+ki)=ires_homo(ki+i01)
+ jres_homo(iistart+ki)=jres_homo(ki+i01)
+ ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+ do k=1,constr_homology
+ odl(k,iistart+ki)=odl(k,ki+i01)
+ sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+ l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+ enddo
+ enddo
+ iistart=iistart+i10-i01
+ endif
+ endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag) then
+ i01=ii
+ liiflag=.true.
+ endif
+ enddo
+ enddo
+ lim_odl=iistart-1
+ endif
+! write (iout,*) "Removing distances completed"
+! call flush(iout)
+ endif ! .not. klapaucjusz
+
+ if (constr_homology.gt.0) call homology_partition
+ write (iout,*) "After homology_partition"
+ call flush(iout)
+ if (constr_homology.gt.0) call init_int_table
+ write (iout,*) "After init_int_table"
+ call flush(iout)
+! endif ! .not. klapaucjusz
+! endif
+! if (constr_homology.gt.0) call homology_partition
+! write (iout,*) "After homology_partition"
+! call flush(iout)
+! if (constr_homology.gt.0) call init_int_table
+! write (iout,*) "After init_int_table"
+! call flush(iout)
+! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+!
+! Print restraints
+!
+#ifdef DEBUG
+ !this debug needs correction
+ if (.not.out_template_restr) return
+!d write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ if(me1.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write (iout,*) "Distance restraints from templates"
+ do ii=1,lim_odl
+ write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
+ ii,ires_homo(ii),jres_homo(ii),&
+ (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
+ ki=1,constr_homology)
+ enddo
+ write (iout,*) "Dihedral angle restraints from templates"
+ do i=nnt+3,nct
+ write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
+ (rad2deg*dih(ki,i),&
+ rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "Virtual-bond angle restraints from templates"
+ do i=nnt+2,nct
+ write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
+ (rad2deg*thetatpl(ki,i),&
+ rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "SC restraints from templates"
+ do i=nnt,nct
+ write(iout,'(i7,100(4f8.2,4x))') i,&
+ (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
+ 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
+ enddo
+ endif
+#endif
+ return
+ end subroutine read_constr_homology
!------------------------------------------------------------------------------
end module io_wham
!-----------------------------------------------------------------------------
projects / unres4.git / blobdiff