X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fio_wham.F90;h=fee1e46c5912644b8b4a9966a5cf0c920860f4fc;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=7ebde6f1ba042e825fdb5e50457b56e761f4e743;hpb=df2469d9ac903d93889867f4e50e9bf6c428c1c6;p=unres4.git diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90 index 7ebde6f..fee1e46 100644 --- a/source/wham/io_wham.F90 +++ b/source/wham/io_wham.F90 @@ -51,6 +51,8 @@ ! 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) @@ -69,6 +71,37 @@ 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') + + #ifndef OLDSCP ! ! 8/9/01 In the newest version SCp interaction constants are read from a file @@ -132,9 +165,9 @@ ! use energy_data use geometry_data, only:nres,deg2rad,c,dc,nres_molec - use control_data, only:iscode + use control_data, only:iscode,dyn_ss 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 ! implicit real*8 (a-h,o-z) @@ -313,6 +346,20 @@ 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 @@ -332,7 +379,7 @@ use geometry_data use energy_data use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor - maxtermd_1,maxtermd_2 !,maxthetyp,maxthetyp1 + maxtermd_1,maxtermd_2,tor_mode,scelemode !,maxthetyp,maxthetyp1 use MD_data !el use MPI_data !el use map_data @@ -365,7 +412,7 @@ character(len=1) :: t1,t2,t3 character(len=1) :: onelett(4) = (/"G","A","P","D"/) character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/) - logical :: lprint + logical :: lprint,SPLIT_FOURIERTOR real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob) character(len=800) :: controlcard character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,& @@ -374,14 +421,19 @@ !el integer ilen !el external ilen character(len=16) :: key - integer :: iparm + 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 + 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 ! @@ -393,10 +445,26 @@ 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) + allocate(ww(max_eneW)) do i=1,n_eneW key = wname(i)(:ilen(wname(i))) @@ -425,11 +493,40 @@ allocate(ww(max_eneW)) wscloc=ww(12) wtor=ww(13) wtor_d=ww(14) + wstrain=ww(15) wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) - wcatcat=ww(44) - wcatprot=ww(43) + 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) +! 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 ! @@ -449,13 +546,34 @@ allocate(ww(max_eneW)) 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(43)=wcatprot - weights(44)=wcatcat + 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 ! el-------- call card_concat(controlcard,.false.) @@ -500,7 +618,7 @@ allocate(ww(max_eneW)) 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)) @@ -534,6 +652,10 @@ allocate(ww(max_eneW)) 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) @@ -580,6 +702,7 @@ allocate(ww(max_eneW)) endif if (.not. allocated(msc)) allocate(msc(ntyp1,5)) if (.not. allocated(restok)) allocate(restok(ntyp1,5)) + if (oldion.eq.1) then do i=1,ntyp_molec(5) read(iion,*) msc(i,5),restok(i,5) @@ -593,6 +716,35 @@ allocate(ww(max_eneW)) 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 !---------------------------------------------------- @@ -738,13 +890,16 @@ allocate(ww(max_eneW)) ! Read the parameters of Utheta determined from ab initio surfaces ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 ! + allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) ! write (iout,*) "tu dochodze" + IF (tor_mode.eq.0) THEN +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! read (ithep,*) nthetyp,ntheterm,ntheterm2,& ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) !---------------------------------------------------- - allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) +! allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) allocate(aa0thet(-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) @@ -864,7 +1019,7 @@ allocate(ww(max_eneW)) ! ! Control printout of the coefficients of virtual-bond-angle potentials ! -do iblock=1,2 + do iblock=1,2 if (lprint) then write (iout,'(//a)') 'Parameter of virtual-bond-angle potential' do i=1,nthetyp+1 @@ -903,8 +1058,92 @@ do iblock=1,2 enddo call flush(iout) endif -enddo + enddo + ELSE +!C here will be the apropriate recalibrating for D-aminoacid + read (ithep,*) nthetyp + allocate(nbend_kcc_Tb(-nthetyp:nthetyp)) + allocate(v1bend_chyb(0:36,-nthetyp:nthetyp)) + do i=-nthetyp+1,nthetyp-1 + read (ithep,*) nbend_kcc_Tb(i) + do j=0,nbend_kcc_Tb(i) + read (ithep,*) ijunk,v1bend_chyb(j,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') & + "Parameters of the valence-only potentials" + do i=-nthetyp+1,nthetyp-1 + write (iout,'(2a)') "Type ",toronelet(i) + do k=0,nbend_kcc_Tb(i) + write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i) + enddo + enddo + endif + 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) @@ -1042,6 +1281,376 @@ enddo 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 + SPLIT_FOURIERTOR = nloctyp.lt.0 + nloctyp = iabs(nloctyp) +#ifdef NEWCORR + if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1)) + print *,"shape",shape(itype2loc) + allocate(iloctyp(-nloctyp:nloctyp)) + allocate(bnew1(3,2,-nloctyp:nloctyp)) + allocate(bnew2(3,2,-nloctyp:nloctyp)) + allocate(ccnew(3,2,-nloctyp:nloctyp)) + allocate(ddnew(3,2,-nloctyp:nloctyp)) + allocate(e0new(3,-nloctyp:nloctyp)) + allocate(eenew(2,2,2,-nloctyp:nloctyp)) + allocate(bnew1tor(3,2,-nloctyp:nloctyp)) + allocate(bnew2tor(3,2,-nloctyp:nloctyp)) + allocate(ccnewtor(3,2,-nloctyp:nloctyp)) + allocate(ddnewtor(3,2,-nloctyp:nloctyp)) + allocate(e0newtor(3,-nloctyp:nloctyp)) + allocate(eenewtor(2,2,2,-nloctyp:nloctyp)) + + read (ifourier,*) (itype2loc(i),i=1,ntyp) + read (ifourier,*) (iloctyp(i),i=0,nloctyp-1) + itype2loc(ntyp1)=nloctyp + iloctyp(nloctyp)=ntyp1 + do i=1,ntyp1 + itype2loc(-i)=-itype2loc(i) + enddo +#else + allocate(iloctyp(-nloctyp:nloctyp)) + allocate(itype2loc(-ntyp1:ntyp1)) + iloctyp(0)=10 + iloctyp(1)=9 + iloctyp(2)=20 + iloctyp(3)=ntyp1 +#endif + do i=1,nloctyp + iloctyp(-i)=-iloctyp(i) + enddo +!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) +!c write (iout,*) "nloctyp",nloctyp, +!c & " iloctyp",(iloctyp(i),i=0,nloctyp) +!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) +!c write (iout,*) "nloctyp",nloctyp, +!c & " iloctyp",(iloctyp(i),i=0,nloctyp) +#ifdef NEWCORR + do i=0,nloctyp-1 +!c write (iout,*) "NEWCORR",i + read (ifourier,*) + do ii=1,3 + do j=1,2 + read (ifourier,*) bnew1(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW1" +!c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2) + do ii=1,3 + do j=1,2 + read (ifourier,*) bnew2(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW2" +!c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*) ccnew(kk,1,i) + read (ifourier,*) ccnew(kk,2,i) + enddo +!c write (iout,*) "NEWCORR CCNEW" +!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*) ddnew(kk,1,i) + read (ifourier,*) ddnew(kk,2,i) + enddo +!c write (iout,*) "NEWCORR DDNEW" +!c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2) + do ii=1,2 + do jj=1,2 + do kk=1,2 + read (ifourier,*) eenew(ii,jj,kk,i) + enddo + enddo + enddo +!c write (iout,*) "NEWCORR EENEW1" +!c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) + do ii=1,3 + read (ifourier,*) e0new(ii,i) + enddo +!c write (iout,*) (e0new(ii,i),ii=1,3) + enddo +!c write (iout,*) "NEWCORR EENEW" + do i=0,nloctyp-1 + do ii=1,3 + ccnew(ii,1,i)=ccnew(ii,1,i)/2 + ccnew(ii,2,i)=ccnew(ii,2,i)/2 + ddnew(ii,1,i)=ddnew(ii,1,i)/2 + ddnew(ii,2,i)=ddnew(ii,2,i)/2 + enddo + enddo + do i=1,nloctyp-1 + do ii=1,3 + bnew1(ii,1,-i)= bnew1(ii,1,i) + bnew1(ii,2,-i)=-bnew1(ii,2,i) + bnew2(ii,1,-i)= bnew2(ii,1,i) + bnew2(ii,2,-i)=-bnew2(ii,2,i) + enddo + do ii=1,3 +!c ccnew(ii,1,i)=ccnew(ii,1,i)/2 +!c ccnew(ii,2,i)=ccnew(ii,2,i)/2 +!c ddnew(ii,1,i)=ddnew(ii,1,i)/2 +!c ddnew(ii,2,i)=ddnew(ii,2,i)/2 + ccnew(ii,1,-i)=ccnew(ii,1,i) + ccnew(ii,2,-i)=-ccnew(ii,2,i) + ddnew(ii,1,-i)=ddnew(ii,1,i) + ddnew(ii,2,-i)=-ddnew(ii,2,i) + enddo + e0new(1,-i)= e0new(1,i) + e0new(2,-i)=-e0new(2,i) + e0new(3,-i)=-e0new(3,i) + do kk=1,2 + eenew(kk,1,1,-i)= eenew(kk,1,1,i) + eenew(kk,1,2,-i)=-eenew(kk,1,2,i) + eenew(kk,2,1,-i)=-eenew(kk,2,1,i) + eenew(kk,2,2,-i)= eenew(kk,2,2,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') "Coefficients of the multibody terms" + do i=-nloctyp+1,nloctyp-1 + write (iout,*) "Type: ",onelet(iloctyp(i)) + write (iout,*) "Coefficients of the expansion of B1" + do j=1,2 + write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of B2" + do j=1,2 + write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of C" + write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3) + write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of D" + write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3) + write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of E" + write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3) + do j=1,2 + do k=1,2 + write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2) + enddo + enddo + enddo + endif + IF (SPLIT_FOURIERTOR) THEN + do i=0,nloctyp-1 +!c write (iout,*) "NEWCORR TOR",i + read (ifourier,*) + do ii=1,3 + do j=1,2 + read (ifourier,*) bnew1tor(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW1 TOR" +!c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2) + do ii=1,3 + do j=1,2 + read (ifourier,*) bnew2tor(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW2 TOR" +!c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*) ccnewtor(kk,1,i) + read (ifourier,*) ccnewtor(kk,2,i) + enddo +!c write (iout,*) "NEWCORR CCNEW TOR" +!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*) ddnewtor(kk,1,i) + read (ifourier,*) ddnewtor(kk,2,i) + enddo +!c write (iout,*) "NEWCORR DDNEW TOR" +!c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2) + do ii=1,2 + do jj=1,2 + do kk=1,2 + read (ifourier,*) eenewtor(ii,jj,kk,i) + enddo + enddo + enddo +!c write (iout,*) "NEWCORR EENEW1 TOR" +!c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) + do ii=1,3 + read (ifourier,*) e0newtor(ii,i) + enddo +!c write (iout,*) (e0newtor(ii,i),ii=1,3) + enddo +!c write (iout,*) "NEWCORR EENEW TOR" + do i=0,nloctyp-1 + do ii=1,3 + ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2 + ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2 + ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2 + ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2 + enddo + enddo + do i=1,nloctyp-1 + do ii=1,3 + bnew1tor(ii,1,-i)= bnew1tor(ii,1,i) + bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i) + bnew2tor(ii,1,-i)= bnew2tor(ii,1,i) + bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i) + enddo + do ii=1,3 + ccnewtor(ii,1,-i)=ccnewtor(ii,1,i) + ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i) + ddnewtor(ii,1,-i)=ddnewtor(ii,1,i) + ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i) + enddo + e0newtor(1,-i)= e0newtor(1,i) + e0newtor(2,-i)=-e0newtor(2,i) + e0newtor(3,-i)=-e0newtor(3,i) + do kk=1,2 + eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i) + eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i) + eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i) + eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') & + "Single-body coefficients of the torsional potentials" + do i=-nloctyp+1,nloctyp-1 + write (iout,*) "Type: ",onelet(iloctyp(i)) + write (iout,*) "Coefficients of the expansion of B1tor" + do j=1,2 + write (iout,'(3hB1(,i1,1h),3f10.5)') & + j,(bnew1tor(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of B2tor" + do j=1,2 + write (iout,'(3hB2(,i1,1h),3f10.5)') & + j,(bnew2tor(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of Ctor" + write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3) + write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of Dtor" + write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3) + write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of Etor" + write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3) + do j=1,2 + do k=1,2 + write (iout,'(1hE,2i1,2f10.5)') & + j,k,(eenewtor(l,j,k,i),l=1,2) + enddo + enddo + enddo + endif + ELSE + do i=-nloctyp+1,nloctyp-1 + do ii=1,3 + do j=1,2 + bnew1tor(ii,j,i)=bnew1(ii,j,i) + enddo + enddo + do ii=1,3 + do j=1,2 + bnew2tor(ii,j,i)=bnew2(ii,j,i) + enddo + enddo + do ii=1,3 + ccnewtor(ii,1,i)=ccnew(ii,1,i) + ccnewtor(ii,2,i)=ccnew(ii,2,i) + ddnewtor(ii,1,i)=ddnew(ii,1,i) + ddnewtor(ii,2,i)=ddnew(ii,2,i) + enddo + enddo + ENDIF !SPLIT_FOURIER_TOR +#else + allocate(ccold(2,2,-nloctyp-1:nloctyp+1)) + allocate(ddold(2,2,-nloctyp-1:nloctyp+1)) + allocate(eeold(2,2,-nloctyp-1:nloctyp+1)) + allocate(b(13,-nloctyp-1:nloctyp+1)) + if (lprint) & + write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)" + do i=0,nloctyp-1 + read (ifourier,*) + read (ifourier,*) (b(ii,i),ii=1,13) + if (lprint) then + write (iout,*) 'Type ',onelet(iloctyp(i)) + write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) + endif + if (i.gt.0) then + b(2,-i)= b(2,i) + b(3,-i)= b(3,i) + b(4,-i)=-b(4,i) + b(5,-i)=-b(5,i) + endif + CCold(1,1,i)= b(7,i) + CCold(2,2,i)=-b(7,i) + CCold(2,1,i)= b(9,i) + CCold(1,2,i)= b(9,i) + CCold(1,1,-i)= b(7,i) + CCold(2,2,-i)=-b(7,i) + CCold(2,1,-i)=-b(9,i) + CCold(1,2,-i)=-b(9,i) + DDold(1,1,i)= b(6,i) + DDold(2,2,i)=-b(6,i) + DDold(2,1,i)= b(8,i) + DDold(1,2,i)= b(8,i) + DDold(1,1,-i)= b(6,i) + DDold(2,2,-i)=-b(6,i) + DDold(2,1,-i)=-b(8,i) + DDold(1,2,-i)=-b(8,i) + EEold(1,1,i)= b(10,i)+b(11,i) + EEold(2,2,i)=-b(10,i)+b(11,i) + EEold(2,1,i)= b(12,i)-b(13,i) + EEold(1,2,i)= b(12,i)+b(13,i) + EEold(1,1,-i)= b(10,i)+b(11,i) + EEold(2,2,-i)=-b(10,i)+b(11,i) + EEold(2,1,-i)=-b(12,i)+b(13,i) + EEold(1,2,-i)=-b(12,i)-b(13,i) + write(iout,*) "TU DOCHODZE" + print *,"JESTEM" + enddo + if (lprint) then + write (iout,*) + write (iout,*) & + "Coefficients of the cumulants (independent of valence angles)" + do i=-nloctyp+1,nloctyp-1 + write (iout,*) 'Type ',onelet(iloctyp(i)) + write (iout,*) 'B1' + write(iout,'(2f10.5)') B(3,i),B(5,i) + write (iout,*) 'B2' + write(iout,'(2f10.5)') B(2,i),B(4,i) + write (iout,*) 'CC' + do j=1,2 + write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i) + enddo + write(iout,*) 'DD' + do j=1,2 + write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i) + enddo + write(iout,*) 'EE' + do j=1,2 + write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) + enddo + enddo + endif +#endif #ifdef CRYST_TOR ! ! Read torsional parameters in old format @@ -1081,6 +1690,8 @@ enddo ! ! Read torsional parameters ! + IF (TOR_MODE.eq.0) THEN + allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) read (itorp,*) ntortyp @@ -1269,6 +1880,118 @@ enddo 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 + read (itorp,*) ntortyp + nkcctyp=ntortyp + write (iout,*) "Valence-torsional parameters read in ntortyp",& + ntortyp + read (itorp,*) (itortyp(i),i=1,ntyp) + write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp) +#ifndef NEWCORR + do i=1,ntyp1 + itype2loc(i)=itortyp(i) + enddo +#endif + do i=-ntyp,-1 + itortyp(i)=-itortyp(-i) + enddo + do i=-ntortyp+1,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 +!C first we read the cos and sin gamma parameters + read (itorp,'(13x,a)') string + write (iout,*) i,j,string + read (itorp,*) & + nterm_kcc(j,i),nterm_kcc_Tb(j,i) +!C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) + do k=1,nterm_kcc(j,i) + do l=1,nterm_kcc_Tb(j,i) + do ll=1,nterm_kcc_Tb(j,i) + read (itorp,*) ii,jj,kk, & + v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) + enddo + enddo + enddo + enddo + enddo + ELSE +#ifdef NEWCORR +!c AL 4/8/16: Calculate coefficients from one-body parameters + ntortyp=nloctyp + allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) + allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1)) + allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1)) + allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1)) + allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1)) + + do i=-ntyp1,ntyp1 + print *,i,itortyp(i) + itortyp(i)=itype2loc(i) + enddo + write (iout,*) & + "Val-tor parameters calculated from cumulant coefficients ntortyp"& + ,ntortyp + do i=-ntortyp+1,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + nterm_kcc(j,i)=2 + nterm_kcc_Tb(j,i)=3 + do k=1,nterm_kcc_Tb(j,i) + do l=1,nterm_kcc_Tb(j,i) + v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)& + +bnew1tor(k,2,i)*bnew2tor(l,2,j) + v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)& + +bnew1tor(k,2,i)*bnew2tor(l,1,j) + enddo + enddo + do k=1,nterm_kcc_Tb(j,i) + do l=1,nterm_kcc_Tb(j,i) +#ifdef CORRCD + v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) & + -ccnewtor(k,2,i)*ddnewtor(l,2,j)) + v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) & + +ccnewtor(k,1,i)*ddnewtor(l,2,j)) +#else + v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) & + -ccnewtor(k,2,i)*ddnewtor(l,2,j)) + v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) & + +ccnewtor(k,1,i)*ddnewtor(l,2,j)) +#endif + enddo + enddo +!c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma) + enddo + enddo +#else + write (iout,*) "TOR_MODE>1 only with NEWCORR" + stop +#endif + ENDIF ! TOR_MODE + if (tor_mode.gt.0 .and. lprint) then +!c Print valence-torsional parameters + write (iout,'(a)') & + "Parameters of the valence-torsional potentials" + do i=-ntortyp+1,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + write (iout,'(3a)') "Type ",toronelet(i),toronelet(j) + write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc" + do k=1,nterm_kcc(j,i) + do l=1,nterm_kcc_Tb(j,i) + do ll=1,nterm_kcc_Tb(j,i) + write (iout,'(3i5,2f15.4)')& + k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) + enddo + enddo + enddo + enddo + enddo + endif + #endif !elwrite(iout,*) "parmread kontrol sc-bb" ! Read of Side-chain backbone correlation parameters @@ -1397,153 +2120,7 @@ enddo enddo enddo endif -! -! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local -! interaction energy of the Gly, Ala, and Pro prototypes. -! - read (ifourier,*) nloctyp -!el write(iout,*)"nloctyp",nloctyp -!el from module energy------- - allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) - allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) - allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) - allocate(cc(2,2,-nloctyp-1:nloctyp+1)) - allocate(dd(2,2,-nloctyp-1:nloctyp+1)) - allocate(ee(2,2,-nloctyp-1:nloctyp+1)) - allocate(ctilde(2,2,-nloctyp-1:nloctyp+1)) - allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor) - do i=1,2 - do ii=-nloctyp-1,nloctyp+1 - b1(i,ii)=0.0d0 - b2(i,ii)=0.0d0 - b1tilde(i,ii)=0.0d0 - do j=1,2 - cc(j,i,ii)=0.0d0 - dd(j,i,ii)=0.0d0 - ee(j,i,ii)=0.0d0 - ctilde(j,i,ii)=0.0d0 - dtilde(j,i,ii)=0.0d0 - enddo - enddo - enddo -!-------------------------------- - allocate(b(13,0:nloctyp)) - - do i=0,nloctyp-1 - read (ifourier,*) - read (ifourier,*) (b(ii,i),ii=1,13) - if (lprint) then - write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) - endif - B1(1,i) = b(3,i) - B1(2,i) = b(5,i) - B1(1,-i) = b(3,i) - B1(2,-i) = -b(5,i) -! b1(1,i)=0.0d0 -! b1(2,i)=0.0d0 - B1tilde(1,i) = b(3,i) - B1tilde(2,i) =-b(5,i) - B1tilde(1,-i) =-b(3,i) - B1tilde(2,-i) =b(5,i) -! b1tilde(1,i)=0.0d0 -! b1tilde(2,i)=0.0d0 - B2(1,i) = b(2,i) - B2(2,i) = b(4,i) - B2(1,-i) =b(2,i) - B2(2,-i) =-b(4,i) - -! b2(1,i)=0.0d0 -! b2(2,i)=0.0d0 - CC(1,1,i)= b(7,i) - CC(2,2,i)=-b(7,i) - CC(2,1,i)= b(9,i) - CC(1,2,i)= b(9,i) - CC(1,1,-i)= b(7,i) - CC(2,2,-i)=-b(7,i) - CC(2,1,-i)=-b(9,i) - CC(1,2,-i)=-b(9,i) -! CC(1,1,i)=0.0d0 -! CC(2,2,i)=0.0d0 -! CC(2,1,i)=0.0d0 -! CC(1,2,i)=0.0d0 - Ctilde(1,1,i)=b(7,i) - Ctilde(1,2,i)=b(9,i) - Ctilde(2,1,i)=-b(9,i) - Ctilde(2,2,i)=b(7,i) - Ctilde(1,1,-i)=b(7,i) - Ctilde(1,2,-i)=-b(9,i) - Ctilde(2,1,-i)=b(9,i) - Ctilde(2,2,-i)=b(7,i) - -! Ctilde(1,1,i)=0.0d0 -! Ctilde(1,2,i)=0.0d0 -! Ctilde(2,1,i)=0.0d0 -! Ctilde(2,2,i)=0.0d0 - DD(1,1,i)= b(6,i) - DD(2,2,i)=-b(6,i) - DD(2,1,i)= b(8,i) - DD(1,2,i)= b(8,i) - DD(1,1,-i)= b(6,i) - DD(2,2,-i)=-b(6,i) - DD(2,1,-i)=-b(8,i) - DD(1,2,-i)=-b(8,i) -! DD(1,1,i)=0.0d0 -! DD(2,2,i)=0.0d0 -! DD(2,1,i)=0.0d0 -! DD(1,2,i)=0.0d0 - Dtilde(1,1,i)=b(6,i) - Dtilde(1,2,i)=b(8,i) - Dtilde(2,1,i)=-b(8,i) - Dtilde(2,2,i)=b(6,i) - Dtilde(1,1,-i)=b(6,i) - Dtilde(1,2,-i)=-b(8,i) - Dtilde(2,1,-i)=b(8,i) - Dtilde(2,2,-i)=b(6,i) - -! Dtilde(1,1,i)=0.0d0 -! Dtilde(1,2,i)=0.0d0 -! Dtilde(2,1,i)=0.0d0 -! Dtilde(2,2,i)=0.0d0 - EE(1,1,i)= b(10,i)+b(11,i) - EE(2,2,i)=-b(10,i)+b(11,i) - EE(2,1,i)= b(12,i)-b(13,i) - EE(1,2,i)= b(12,i)+b(13,i) - EE(1,1,-i)= b(10,i)+b(11,i) - EE(2,2,-i)=-b(10,i)+b(11,i) - EE(2,1,-i)=-b(12,i)+b(13,i) - EE(1,2,-i)=-b(12,i)-b(13,i) - -! ee(1,1,i)=1.0d0 -! ee(2,2,i)=1.0d0 -! ee(2,1,i)=0.0d0 -! ee(1,2,i)=0.0d0 -! ee(2,1,i)=ee(1,2,i) - enddo - if (lprint) then - do i=1,nloctyp - write (iout,*) 'Type',i - write (iout,*) 'B1' -! write (iout,'(f10.5)') B1(:,i) - write(iout,*) B1(1,i),B1(2,i) - write (iout,*) 'B2' -! write (iout,'(f10.5)') B2(:,i) - write(iout,*) B2(1,i),B2(2,i) - write (iout,*) 'CC' - do j=1,2 - write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i) - enddo - write(iout,*) 'DD' - do j=1,2 - write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i) - enddo - write(iout,*) 'EE' - do j=1,2 - write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i) - enddo - enddo - endif ! ! Read electrostatic-interaction parameters ! @@ -1632,6 +2209,7 @@ enddo !---------------------- GB or BP potential ----------------------------- case (3:4) ! 30 do i=1,ntyp + if (scelemode.eq.0) then do i=1,ntyp read (isidep,*)(eps(i,j),j=i,ntyp) enddo @@ -1658,6 +2236,86 @@ enddo write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),& chip(i),alp(i),i=1,ntyp) endif + else + allocate(icharge(ntyp1)) +! print *,ntyp,icharge(i) + icharge(:)=0 + read (isidep,*) (icharge(i),i=1,ntyp) + print *,ntyp,icharge(i) +! if(.not.allocated(eps)) allocate(eps(-ntyp + write (2,*) "icharge",(icharge(i),i=1,ntyp) + allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp)) + allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp)) + allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp)) + allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp)) + allocate(epsintab(ntyp,ntyp)) + allocate(dtail(2,ntyp,ntyp)) + print *,"control line 1" + allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp)) + allocate(wqdip(2,ntyp,ntyp)) + allocate(wstate(4,ntyp,ntyp)) + allocate(dhead(2,2,ntyp,ntyp)) + allocate(nstate(ntyp,ntyp)) + allocate(debaykap(ntyp,ntyp)) + print *,"control line 2" + if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) + if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp) + + do i=1,ntyp + do j=1,i +! write (*,*) "Im in ALAB", i, " ", j + read(isidep,*) & + eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & + (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & + chis(i,j),chis(j,i), & + nstate(i,j),(wstate(k,i,j),k=1,4), & + dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& + dtail(1,i,j),dtail(2,i,j), & + epshead(i,j),sig0head(i,j), & + rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & + alphapol(i,j),alphapol(j,i), & + (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) +! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) + END DO + END DO + DO i = 1, ntyp + DO j = i+1, ntyp + eps(i,j) = eps(j,i) + sigma(i,j) = sigma(j,i) + sigmap1(i,j)=sigmap1(j,i) + sigmap2(i,j)=sigmap2(j,i) + sigiso1(i,j)=sigiso1(j,i) + sigiso2(i,j)=sigiso2(j,i) +! print *,"ATU",sigma(j,i),sigma(i,j),i,j + nstate(i,j) = nstate(j,i) + dtail(1,i,j) = dtail(1,j,i) + dtail(2,i,j) = dtail(2,j,i) + DO k = 1, 4 + alphasur(k,i,j) = alphasur(k,j,i) + wstate(k,i,j) = wstate(k,j,i) + alphiso(k,i,j) = alphiso(k,j,i) + END DO + + dhead(2,1,i,j) = dhead(1,1,j,i) + dhead(2,2,i,j) = dhead(1,2,j,i) + dhead(1,1,i,j) = dhead(2,1,j,i) + dhead(1,2,i,j) = dhead(2,2,j,i) + + epshead(i,j) = epshead(j,i) + sig0head(i,j) = sig0head(j,i) + + DO k = 1, 2 + wqdip(k,i,j) = wqdip(k,j,i) + END DO + + wquad(i,j) = wquad(j,i) + epsintab(i,j) = epsintab(j,i) + debaykap(i,j)=debaykap(j,i) +! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j) + END DO + END DO + endif + ! goto 50 !--------------------- GBV potential ----------------------------------- case (5) @@ -1688,18 +2346,26 @@ enddo !el from module energy - COMMON.INTERACT------- ! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) - allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) + allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) - allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) + 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 bb_aq(i,j)=0.0D0 aa_lip(i,j)=0.0d0 bb_lip(i,j)=0.0d0 + if (scelemode.eq.0) then chi(i,j)=0.0D0 sigma(i,j)=0.0D0 r0(i,j)=0.0D0 + endif enddo enddo !-------------------------------- @@ -1710,6 +2376,7 @@ enddo epslip(i,j)=epslip(j,i) enddo enddo + if (scelemode.eq.0) then do i=1,ntyp do j=i,ntyp sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2) @@ -1718,6 +2385,7 @@ enddo rs0(j,i)=rs0(i,j) enddo enddo + endif if (lprint) write (iout,'(/a/10x,7a/72(1h-))') & 'Working parameters of the SC interactions:',& ' a ',' b ',' augm ',' sigma ',' r0 ',& @@ -1727,6 +2395,7 @@ enddo epsij=eps(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) + print *,"rrij",rrij else rrij=rr0(i)+rr0(j) endif @@ -1747,7 +2416,7 @@ enddo bb_lip(i,j)=-sigeps*epsijlip*rrij aa_lip(j,i)=aa_lip(i,j) bb_lip(j,i)=bb_lip(i,j) - if (ipot.gt.2) then + if ((ipot.gt.2).and. (scelemode.eq.0))then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 sigii1=sigii(i) @@ -1784,6 +2453,271 @@ enddo 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 @@ -1791,6 +2725,126 @@ enddo 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))) +! 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=1,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) + 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) + 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 + #ifdef CLUSTER ! ! Define the SC-p interaction constants @@ -1802,6 +2856,48 @@ enddo 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" ! @@ -1852,10 +2948,24 @@ enddo 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. ! @@ -1866,14 +2976,28 @@ enddo ! 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 @@ -1952,9 +3076,12 @@ enddo !----------------------------------------------------------------------------- subroutine read_general_data(*) - use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions - use energy_data, only:distchainmax - use geometry_data, only:boxxsize,boxysize,boxzsize + use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,& + scelemode,TUBEmode,tor_mode,energy_dec + + use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss + use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,& + bordtubebot,tubebufthick,buftubebot,buftubetop ! implicit none ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" @@ -2034,8 +3161,30 @@ enddo call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) + 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) + + if (TUBEmode.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + 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 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0 @@ -2656,6 +3805,7 @@ enddo 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) @@ -2758,8 +3908,9 @@ enddo !-------------------------------------------------------------------------------- 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' @@ -2770,10 +3921,11 @@ enddo ! 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 @@ -2788,16 +3940,26 @@ enddo mnum=molnum(i) iti=itype(i,mnum) if (iti.eq.ntyp1) then + if (itype(i-1,molnum(i-1)).eq.ntyp1) 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)