X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.F90;h=3dc6fa1fb4a2b3fbc17f17c97b8405907dfcec0c;hb=3d3780d89d2885ee8700324ef832ef07ee35c45b;hp=87ddf94beea67828febd369586e764f506f5e616;hpb=df2469d9ac903d93889867f4e50e9bf6c428c1c6;p=unres4.git diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90 index 87ddf94..3dc6fa1 100644 --- a/source/unres/io_config.F90 +++ b/source/unres/io_config.F90 @@ -734,19 +734,21 @@ 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,LaTeX + logical :: lprint,LaTeX,SPLIT_FOURIERTOR real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob) real(kind=8),dimension(13) :: buse character(len=3) :: lancuch !,ucase !el local variables - integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm - integer :: maxinter,junk,kk,ii,ncatprotparm + integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj + integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,& dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,& res1,epsijlip,epspeptube,epssctube,sigmapeptube, & sigmasctube - integer :: ichir1,ichir2 + integer :: ichir1,ichir2,ijunk + character*3 string + ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) @@ -865,8 +867,13 @@ enddo enddo endif + + + + if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5))) + if (oldion.eq.1) then do i=1,ntyp_molec(5) - read(iion,*) msc(i,5),restok(i,5) + read(iion,*) msc(i,5),restok(i,5),ichargecat(i) print *,msc(i,5),restok(i,5) enddo ip(5)=0.2 @@ -877,6 +884,23 @@ 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 (iion,*) (vcatprm(k),k=1,ncatprotpram) !---------------------------------------------------- allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp)) @@ -1017,6 +1041,7 @@ ! Read the parameters of Utheta determined from ab initio surfaces ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 ! + IF (tor_mode.eq.0) THEN read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,& ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) @@ -1209,6 +1234,29 @@ enddo call flush(iout) endif + ELSE +!C here will be the apropriate recalibrating for D-aminoacid + read (ithep,*,end=121,err=121) nthetyp + allocate(nbend_kcc_Tb(-nthetyp:nthetyp)) + allocate(v1bend_chyb(0:36,-nthetyp:nthetyp)) + do i=-nthetyp+1,nthetyp-1 + read (ithep,*,end=121,err=121) nbend_kcc_Tb(i) + do j=0,nbend_kcc_Tb(i) + read (ithep,*,end=121,err=121) 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 + write (2,*) "Start reading THETA_PDB",ithep_pdb do i=1,ntyp ! write (2,*) 'i=',i @@ -1501,6 +1549,444 @@ #endif close(irotam) + +!C +!C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local +!C interaction energy of the Gly, Ala, and Pro prototypes. +!C + read (ifourier,*) nloctyp + SPLIT_FOURIERTOR = nloctyp.lt.0 + nloctyp = iabs(nloctyp) +!C allocate(b1(2,nres)) !(2,-maxtor:maxtor) +!C allocate(b2(2,nres)) !(2,-maxtor:maxtor) +!C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor) +!C allocate(ctilde(2,2,nres)) +!C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor) +!C allocate(gtb1(2,nres)) +!C allocate(gtb2(2,nres)) +!C allocate(cc(2,2,nres)) +!C allocate(dd(2,2,nres)) +!C allocate(ee(2,2,nres)) +!C allocate(gtcc(2,2,nres)) +!C allocate(gtdd(2,2,nres)) +!C allocate(gtee(2,2,nres)) + +#ifdef NEWCORR + allocate(itype2loc(-ntyp1:ntyp1)) + 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,*,end=115,err=115) (itype2loc(i),i=1,ntyp) + read (ifourier,*,end=115,err=115) (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,*,end=115,err=115) + do ii=1,3 + do j=1,2 + read (ifourier,*,end=115,err=115) 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,*,end=115,err=115) 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,*,end=115,err=115) ccnew(kk,1,i) + read (ifourier,*,end=115,err=115) 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,*,end=115,err=115) ddnew(kk,1,i) + read (ifourier,*,end=115,err=115) 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,*,end=115,err=115) 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,*,end=115,err=115) 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,*,end=115,err=115) + do ii=1,3 + do j=1,2 + read (ifourier,*,end=115,err=115) 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,*,end=115,err=115) 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,*,end=115,err=115) ccnewtor(kk,1,i) + read (ifourier,*,end=115,err=115) 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,*,end=115,err=115) ddnewtor(kk,1,i) + read (ifourier,*,end=115,err=115) 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,*,end=115,err=115) 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,*,end=115,err=115) 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,*,end=115,err=115) + read (ifourier,*,end=115,err=115) (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 +!c B1(1,i) = b(3) +!c B1(2,i) = b(5) +!c B1(1,-i) = b(3) +!c B1(2,-i) = -b(5) +!c b1(1,i)=0.0d0 +!c b1(2,i)=0.0d0 +!c B1tilde(1,i) = b(3) +!c! B1tilde(2,i) =-b(5) +!c! B1tilde(1,-i) =-b(3) +!c! B1tilde(2,-i) =b(5) +!c! b1tilde(1,i)=0.0d0 +!c b1tilde(2,i)=0.0d0 +!c B2(1,i) = b(2) +!c B2(2,i) = b(4) +!c B2(1,-i) =b(2) +!c B2(2,-i) =-b(4) +!cc B1tilde(1,i) = b(3,i) +!cc B1tilde(2,i) =-b(5,i) +!c B1tilde(1,-i) =-b(3,i) +!c B1tilde(2,-i) =b(5,i) +!cc b1tilde(1,i)=0.0d0 +!cc b1tilde(2,i)=0.0d0 +!cc B2(1,i) = b(2,i) +!cc B2(2,i) = b(4,i) +!c B2(1,-i) =b(2,i) +!c B2(2,-i) =-b(4,i) + +!c b2(1,i)=0.0d0 +!c b2(2,i)=0.0d0 + 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) +!c CC(1,1,i)=0.0d0 +!c CC(2,2,i)=0.0d0 +!c CC(2,1,i)=0.0d0 +!c CC(1,2,i)=0.0d0 +!c Ctilde(1,1,i)= CCold(1,1,i) +!c Ctilde(1,2,i)= CCold(1,2,i) +!c Ctilde(2,1,i)=-CCold(2,1,i) +!c Ctilde(2,2,i)=-CCold(2,2,i) +!c CC(1,1,i)=0.0d0 +!c CC(2,2,i)=0.0d0 +!c CC(2,1,i)=0.0d0 +!c CC(1,2,i)=0.0d0 +!c Ctilde(1,1,i)= CCold(1,1,i) +!c Ctilde(1,2,i)= CCold(1,2,i) +!c Ctilde(2,1,i)=-CCold(2,1,i) +!c Ctilde(2,2,i)=-CCold(2,2,i) + +!c Ctilde(1,1,i)=0.0d0 +!c Ctilde(1,2,i)=0.0d0 +!c Ctilde(2,1,i)=0.0d0 +!c Ctilde(2,2,i)=0.0d0 + 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) +!c DD(1,1,i)=0.0d0 +!c DD(2,2,i)=0.0d0 +!c DD(2,1,i)=0.0d0 +!c DD(1,2,i)=0.0d0 +!c Dtilde(1,1,i)= DD(1,1,i) +!c Dtilde(1,2,i)= DD(1,2,i) +!c Dtilde(2,1,i)=-DD(2,1,i) +!c Dtilde(2,2,i)=-DD(2,2,i) + +!c Dtilde(1,1,i)=0.0d0 +!c Dtilde(1,2,i)=0.0d0 +!c Dtilde(2,1,i)=0.0d0 +!c Dtilde(2,2,i)=0.0d0 + 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" +!c ee(1,1,i)=1.0d0 +!c ee(2,2,i)=1.0d0 +!c ee(2,1,i)=0.0d0 +!c ee(1,2,i)=0.0d0 +!c ee(2,1,i)=ee(1,2,i) + 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 @@ -1537,6 +2023,7 @@ ! ! Read torsional parameters ! + IF (TOR_MODE.eq.0) THEN allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) read (itorp,*,end=113,err=113) ntortyp !el from energy module--------- @@ -1722,6 +2209,119 @@ 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,*,end=121,err=121) ntortyp + nkcctyp=ntortyp + write (iout,*) "Valence-torsional parameters read in ntortyp",& + ntortyp + read (itorp,*,end=121,err=121) (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)',end=121,err=121) string + write (iout,*) i,j,string + read (itorp,*,end=121,err=121) & + 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,*,end=121,err=121) 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 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) read (itorp_nucl,*,end=113,err=113) ntortyp_nucl @@ -1930,139 +2530,6 @@ ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local ! interaction energy of the Gly, Ala, and Pro prototypes. ! - - if (lprint) then - write (iout,*) - write (iout,*) "Coefficients of the cumulants" - endif - read (ifourier,*) nloctyp -!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) -! el - b1(1,:)=0.0d0 - b1(2,:)=0.0d0 -!-------------------------------- - - do i=0,nloctyp-1 - read (ifourier,*,end=115,err=115) - read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13) - if (lprint) then - write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13) - endif - B1(1,i) = buse(3) - B1(2,i) = buse(5) - B1(1,-i) = buse(3) - B1(2,-i) = -buse(5) -! buse1(1,i)=0.0d0 -! buse1(2,i)=0.0d0 - B1tilde(1,i) = buse(3) - B1tilde(2,i) =-buse(5) - B1tilde(1,-i) =-buse(3) - B1tilde(2,-i) =buse(5) -! buse1tilde(1,i)=0.0d0 -! buse1tilde(2,i)=0.0d0 - B2(1,i) = buse(2) - B2(2,i) = buse(4) - B2(1,-i) =buse(2) - B2(2,-i) =-buse(4) - -! buse2(1,i)=0.0d0 -! buse2(2,i)=0.0d0 - CC(1,1,i)= buse(7) - CC(2,2,i)=-buse(7) - CC(2,1,i)= buse(9) - CC(1,2,i)= buse(9) - CC(1,1,-i)= buse(7) - CC(2,2,-i)=-buse(7) - CC(2,1,-i)=-buse(9) - CC(1,2,-i)=-buse(9) -! 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)=buse(7) - Ctilde(1,2,i)=buse(9) - Ctilde(2,1,i)=-buse(9) - Ctilde(2,2,i)=buse(7) - Ctilde(1,1,-i)=buse(7) - Ctilde(1,2,-i)=-buse(9) - Ctilde(2,1,-i)=buse(9) - Ctilde(2,2,-i)=buse(7) - -! 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)= buse(6) - DD(2,2,i)=-buse(6) - DD(2,1,i)= buse(8) - DD(1,2,i)= buse(8) - DD(1,1,-i)= buse(6) - DD(2,2,-i)=-buse(6) - DD(2,1,-i)=-buse(8) - DD(1,2,-i)=-buse(8) -! 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)=buse(6) - Dtilde(1,2,i)=buse(8) - Dtilde(2,1,i)=-buse(8) - Dtilde(2,2,i)=buse(6) - Dtilde(1,1,-i)=buse(6) - Dtilde(1,2,-i)=-buse(8) - Dtilde(2,1,-i)=buse(8) - Dtilde(2,2,-i)=buse(6) - -! 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)= buse(10)+buse(11) - EE(2,2,i)=-buse(10)+buse(11) - EE(2,1,i)= buse(12)-buse(13) - EE(1,2,i)= buse(12)+buse(13) - EE(1,1,-i)= buse(10)+buse(11) - EE(2,2,-i)=-buse(10)+buse(11) - EE(2,1,-i)=-buse(12)+buse(13) - EE(1,2,-i)=-buse(12)-buse(13) - -! 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,*) B1(1,i),B1(2,i) - write (iout,*) 'B2' - 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 ! @@ -2207,23 +2674,45 @@ allocate(wstate(4,ntyp,ntyp)) allocate(dhead(2,2,ntyp,ntyp)) allocate(nstate(ntyp,ntyp)) + allocate(debaykap(ntyp,ntyp)) + 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) -! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) + eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii + (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii + chis(i,j),chis(j,i), & !2 w tej linii + nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy + dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii + dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini + epshead(i,j),sig0head(i,j), & ! 2 w tej linii + rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii + alphapol(i,j),alphapol(j,i), & ! 2 w tej linii + (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii + IF ((LaTeX).and.(i.gt.24)) then + write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), & + eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii + (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii + chis(i,j),chis(j,i) !2 w tej linii + endif + print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j) + END DO + END DO + do i=1,ntyp + do j=1,i + IF ((LaTeX).and.(i.gt.24)) then + write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), & + dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii + dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini + epshead(i,j),sig0head(i,j), & ! 2 w tej linii + rborn(i,j),rborn(j,i), & ! 3 w tej linii + alphapol(i,j),alphapol(j,i), & ! 2 w tej linii + (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii + endif END DO END DO DO i = 1, ntyp @@ -2258,10 +2747,13 @@ 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) @@ -2353,7 +2845,7 @@ sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_aq(i,j)=epsij*rrij*rrij - print *,"ADASKO",epsij,rrij,expon +! print *,"ADASKO",epsij,rrij,expon bb_aq(i,j)=-sigeps*epsij*rrij aa_aq(j,i)=aa_aq(i,j) bb_aq(j,i)=bb_aq(i,j) @@ -2762,6 +3254,158 @@ ! v1ss=0.0d0 ! v2ss=0.0d0 ! v3ss=0.0d0 + +! Ions by Aga + + allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(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) + + + if (.not.allocated(ichargecat)) 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),alphapolcat(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) +! if (i.eq.1) then +! write (iout,*) 'i= ', i, ' j= ', j +! write (iout,*) 'epsi0= ', epscat(1,j) +! write (iout,*) 'sigma0= ', sigmacat(1,j) +! write (iout,*) 'chi(i,j)= ', chicat(1,j) +! write (iout,*) 'chi(j,i)= ', chicat(j,1) +! write (iout,*) 'chip(i,j)= ', chippcat(1,j) +! write (iout,*) 'chip(j,i)= ', chippcat(j,1) +! 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,*) 'chis(i,j)= ', chiscat(1,j) +! write (iout,*) 'chis(j,i)= ', chiscat(j,1) +! write (iout,*) 'dhead= ', dheadcat(1,1,1,j) +! write (iout,*) 'a1= ', rborncat(j,1) +! write (iout,*) 'a2= ', rborncat(1,j) +! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1) +! write (iout,*) 'alphapol1= ', alphapolcat(1,j) +! write (iout,*) 'w1= ', wqdipcat(1,1,j) +! write (iout,*) 'w2= ', wqdipcat(2,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 + + 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 + if(me.eq.king) then write (iout,'(/a)') "Disulfide bridge parameters:" @@ -2815,6 +3459,8 @@ 118 write (iout,*) "Error reading SCp interaction parameters." goto 999 119 write (iout,*) "Error reading SCCOR parameters" + go to 999 + 121 write (iout,*) "Error in Czybyshev parameters" 999 continue #ifdef MPI call MPI_Finalize(Ierror) @@ -2876,7 +3522,7 @@ logical :: lprn=.true.,fail real(kind=8),dimension(3) :: e1,e2,e3 real(kind=8) :: dcj,efree_temp - character(len=3) :: seq,res + character(len=3) :: seq,res,res2 character(len=5) :: atom character(len=80) :: card real(kind=8),dimension(3,20) :: sccor @@ -2968,6 +3614,8 @@ if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp ! Fish out the ATOM cards. ! write(iout,*) 'card',card(1:20) +! print *,"ATU ",card(1:6), CARD(3:6) +! print *,card if (index(card(1:4),'ATOM').gt.0) then read (card(12:16),*) atom ! write (iout,*) "! ",atom," !",ires @@ -2989,9 +3637,9 @@ enddo else call sccenter(ires_old,iii,sccor) - endif + endif !unres_pdb iii=0 - endif + endif !ind_pdb ! Start new residue. if (res.eq.'Cl-' .or. res.eq.'Na+') then ires=ires_old @@ -3003,7 +3651,7 @@ ishift=ishift-1 itype(1,1)=ntyp1 nres_molec(molecule)=nres_molec(molecule)+1 - endif + endif ! Gly ires=ires-ishift+ishift1 ires_old=ires ! write (iout,*) "ishift",ishift," ires",ires,& @@ -3022,7 +3670,7 @@ ishift=ishift-(ires-ishift+ishift1-ires_old-1) ires=ires-ishift+ishift1 ires_old=ires - endif + endif ! Na Cl ! print *,'atom',ires,atom if (res.eq.'ACE' .or. res.eq.'NHE') then itype(ires,1)=10 @@ -3034,7 +3682,8 @@ ! nres_molec(molecule)=nres_molec(molecule)+1 else molecule=2 - itype(ires,molecule)=rescode(ires,res(2:3),0,molecule) + res2=res(2:3) + itype(ires,molecule)=rescode(ires,res2,0,molecule) ! nres_molec(molecule)=nres_molec(molecule)+1 read (card(19:19),'(a1)') sugar isugar=sugarcode(sugar,ires) @@ -3045,11 +3694,11 @@ ! print *,"ires=",ires,istype(ires) ! endif - endif - endif + endif ! atom.eq.CA + endif !ACE else ires=ires-ishift+ishift1 - endif + endif !ires_old ! write (iout,*) "ires_old",ires_old," ires",ires if (card(27:27).eq."A" .or. card(27:27).eq."B") then ! ishift1=ishift1+1 @@ -3134,6 +3783,7 @@ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) endif endif +! print *,"IONS",ions,card(1:6) else if ((ions).and.(card(1:6).eq.'HETATM')) then if (firstion.eq.0) then firstion=1 @@ -3143,10 +3793,10 @@ enddo else call sccenter(ires,iii,sccor) - endif - endif + endif ! unres_pdb + endif !firstion read (card(12:16),*) atom - print *,"HETATOM", atom +! print *,"HETATOM", atom read (card(18:20),'(a3)') res if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.& (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') & @@ -3160,7 +3810,7 @@ res=res(2:3)//' ' itype(ires,molecule)=rescode(ires,res,0,molecule) read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) - endif + endif! NA endif !atom enddo 10 write (iout,'(a,i5)') ' Number of residues found: ',ires @@ -3234,7 +3884,8 @@ ! e2(3)=0.0d0 ! endif do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) +! c(j,i)=c(j,i+1)-1.9d0*e2(j) + c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo ! else !unres_pdb do j=1,3 @@ -3349,7 +4000,8 @@ e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)-1.9d0*e2(j) +! c(j,1)=c(j,2)-1.9d0*e2(j) + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0) enddo else do j=1,3 @@ -3442,6 +4094,18 @@ istype(i)=istype_temp(i) enddo enddo + if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then +! I have only ions now dummy atoms in the system + molnum(1)=5 + itype(1,5)=itype(2,5) + itype(1,1)=0 + do i=2,nres + itype(i,5)=itype(i+1,5) + enddo + itype(nres,5)=0 + nres=nres-1 + nres_molec(1)=nres_molec(1)-1 + endif ! if (itype(1,1).eq.ntyp1) then ! nsup=nsup-1 ! nstart_sup=2 @@ -3477,7 +4141,7 @@ enddo endif -! print *,seqalingbegin,nres + print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -3513,13 +4177,17 @@ allocate(dc_norm(3,0:2*nres+2)) dc_norm(:,:)=0.d0 endif - + write(iout,*) "before int_from_cart" call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) + write(iout,*) "after int_from_cart" + + do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo + write(iout,*) "after thetaref" ! do i=1,2*nres ! vbld_inv(i)=0.d0 ! vbld(i)=0.d0 @@ -3549,6 +4217,13 @@ ! enddo ! enddo ! +! do i=1,2*nres +! do j=1,3 +! chomo(j,i,k)=c(j,i) +! enddo +! enddo +! write(iout,*) "after chomo" + if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm) if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym) if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym) @@ -3712,7 +4387,8 @@ real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,& integer i - + usampl=.false. ! the default value of usample should be 0 +! write(iout,*) "KURWA2",usampl nglob_csa=0 eglob_csa=1d99 nmin_csa=0 @@ -3734,6 +4410,10 @@ call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0) call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0) call reada(controlcard,'DRMS',drms,0.1D0) + call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) + 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(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 @@ -3758,10 +4438,15 @@ call readi(controlcard,'SHIELD',shield_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "shield_mode",shield_mode + 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 + !C Varibles set size of box with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 protein=index(controlcard,"PROTEIN").gt.0 ions=index(controlcard,"IONS").gt.0 + call readi(controlcard,'OLDION',oldion,1) nucleic=index(controlcard,"NUCLEIC").gt.0 write (iout,*) "with_theta_constr ",with_theta_constr AFMlog=(index(controlcard,'AFM')) @@ -3779,6 +4464,8 @@ ! elemode = 0 is orignal UNRES electrostatics ! elemode = 1 is "Momo" potentials in progress ! elemode = 2 is in development EVALD + + write (iout,*) TUBEmode,"TUBEMODE" if (TUBEmode.gt.0) then call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) @@ -3793,7 +4480,7 @@ endif ! CUTOFFF ON ELECTROSTATICS - 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 ! Lipidic parameters @@ -3934,6 +4621,7 @@ ! enddo buff_shield=1.0d0 endif + itime_mat=0 return end subroutine read_control !----------------------------------------------------------------------------- @@ -4078,6 +4766,7 @@ rest = index(controlcard,"REST").gt.0 tbf = index(controlcard,"TBF").gt.0 usampl = index(controlcard,"USAMPL").gt.0 +! write(iout,*) "KURWA",usampl mdpdb = index(controlcard,"MDPDB").gt.0 call reada(controlcard,"T_BATH",t_bath,300.0d0) call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) @@ -4684,6 +5373,8 @@ open (itube,file=tubename,status='old') call getenv_loc('IONPAR',ionname) open (iion,file=ionname,status='old') + call getenv_loc('IONPAR_NUCL',ionnuclname) + open (iionnucl,file=ionnuclname,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -4743,6 +5434,8 @@ open (itube,file=tubename,status='old',action='read') call getenv_loc('IONPAR',ionname) open (iion,file=ionname,status='old',action='read') + call getenv_loc('IONPAR_NUCL',ionnuclname) + open (iionnucl,file=ionnuclname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) @@ -5345,6 +6038,332 @@ return end subroutine write_stat_thread !----------------------------------------------------------------------------- + subroutine readpdb_template(k) +! Read the PDB file for read_constr_homology with read2sigma +! and convert the peptide geometry into virtual-chain geometry. +! implicit none +! include 'DIMENSIONS' +! include 'COMMON.LOCAL' +! include 'COMMON.VAR' +! include 'COMMON.CHAIN' +! include 'COMMON.INTERACT' +! include 'COMMON.IOUNITS' +! include 'COMMON.GEO' +! include 'COMMON.NAMES' +! include 'COMMON.CONTROL' +! include 'COMMON.FRAG' +! include 'COMMON.SETUP' + use compare_data, only:nhfrag,nbfrag + integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, & + ishift_pdb,ires_ca + logical lprn /.false./,fail + real(kind=8), dimension (3):: e1,e2,e3 + real(kind=8) :: dcj,efree_temp + character*3 seq,res + character*5 atom + character*80 card + real(kind=8), dimension (3,20) :: sccor +! integer rescode + integer, dimension (:), allocatable :: iterter + if(.not.allocated(iterter))allocate(iterter(nres)) + do i=1,nres + iterter(i)=0 + enddo + ibeg=1 + ishift1=0 + ishift=0 + write (2,*) "UNRES_PDB",unres_pdb + ires=0 + ires_old=0 + iii=0 + lsecondary=.false. + nhfrag=0 + nbfrag=0 + do + read (ipdbin,'(a80)',end=10) card + if (card(:3).eq.'END') then + goto 10 + else if (card(:3).eq.'TER') then +! End current chain + ires_old=ires+2 + itype(ires_old-1,1)=ntyp1 + iterter(ires_old-1)=1 + itype(ires_old,1)=ntyp1 + iterter(ires_old)=1 + ibeg=2 +! write (iout,*) "Chain ended",ires,ishift,ires_old + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + endif +! Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + read (card(12:16),*) atom +! write (iout,*) "! ",atom," !",ires +! if (atom.eq.'CA' .or. atom.eq.'CH3') then + read (card(23:26),*) ires + read (card(18:20),'(a3)') res +! write (iout,*) "ires",ires,ires-ishift+ishift1, +! & " ires_old",ires_old +! write (iout,*) "ishift",ishift," ishift1",ishift1 +! write (iout,*) "IRES",ires-ishift+ishift1,ires_old + if (ires-ishift+ishift1.ne.ires_old) then +! Calculate the CM of the preceding residue. + if (ibeg.eq.0) then + if (unres_pdb) then + do j=1,3 + dc(j,ires_old)=sccor(j,iii) + enddo + else + call sccenter(ires_old,iii,sccor) + endif + iii=0 + endif +! Start new residue. + if (res.eq.'Cl-' .or. res.eq.'Na+') then + ires=ires_old + cycle + else if (ibeg.eq.1) then +! write (iout,*) "BEG ires",ires + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1,1)=ntyp1 + endif + ires=ires-ishift+ishift1 + ires_old=ires +! write (iout,*) "ishift",ishift," ires",ires, +! & " ires_old",ires_old +! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift + ibeg=0 + else if (ibeg.eq.2) then +! Start a new chain + ishift=-ires_old+ires-1 + ires=ires_old+1 +! write (iout,*) "New chain started",ires,ishift + ibeg=0 + else + ishift=ishift-(ires-ishift+ishift1-ires_old-1) + ires=ires-ishift+ishift1 + ires_old=ires + endif + if (res.eq.'ACE' .or. res.eq.'NHE') then + itype(ires,1)=10 + else + itype(ires,1)=rescode(ires,res,0,1) + endif + else + ires=ires-ishift+ishift1 + endif +! write (iout,*) "ires_old",ires_old," ires",ires +! if (card(27:27).eq."A" .or. card(27:27).eq."B") then +! ishift1=ishift1+1 +! endif +! write (2,*) "ires",ires," res ",res," ity",ity + if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & + res.eq.'NHE'.and.atom(:2).eq.'HN') then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3) +#ifdef DEBUG + write (iout,'(2i3,2x,a,3f8.3)') & + ires,itype(ires,1),res,(c(j,ires),j=1,3) +#endif + iii=iii+1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo + if (ishift.ne.0) then + ires_ca=ires+ishift-ishift1 + else + ires_ca=ires + endif +! write (*,*) card(23:27),ires,itype(ires) + else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.& + atom.ne.'N' .and. atom.ne.'C' .and.& + atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.& + atom.ne.'OXT' .and. atom(:2).ne.'3H') then +! write (iout,*) "sidechain ",atom + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + endif + enddo + 10 if(me.eq.king.or..not.out1file) & + write (iout,'(a,i5)') ' Nres: ',ires +! Calculate dummy residue coordinates inside the "chain" of a multichain +! system + nres=ires + do i=2,nres-1 +! write (iout,*) i,itype(i),itype(i+1) + if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then + if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then +! 16/01/2014 by Adasko: Adding to dummy atoms in the chain +! first is connected prevous chain (itype(i+1).eq.ntyp1)=true +! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false + if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(i-3,i-2,i-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif !fail + do j=1,3 + c(j,i)=c(j,i-1)-1.9d0*e2(j) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + enddo + endif !unres_pdb + else !itype(i+1).eq.ntyp1 + if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(i+1,i+2,i+3,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,i)=c(j,i+1)-1.9d0*e2(j) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + enddo + endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 + enddo +! Calculate the CM of the last side chain. + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + nsup=nres + nstart_sup=1 + if (itype(nres,1).ne.10) then + nres=nres+1 + itype(nres,1)=ntyp1 + if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,nres)=c(j,nres-1)-1.9d0*e2(j) + enddo + else + do j=1,3 + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1,1).eq.ntyp1) then + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)-1.9d0*e2(j) + enddo + else + do j=1,3 + dcj=(c(j,4)-c(j,3))/2.0 + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif + endif +! Copy the coordinates to reference coordinates +! do i=1,2*nres +! do j=1,3 +! cref(j,i)=c(j,i) +! enddo +! enddo +! Calculate internal coordinates. + if (out_template_coord) then + write (iout,'(/a)') & + "Cartesian coordinates of the reference structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & + "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')& + restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),& + (c(j,ires+nres),j=1,3) + enddo + endif +! Calculate internal coordinates. + call int_from_cart(.true.,out_template_coord) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo +! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), +! & vbld_inv(i+nres) + enddo + do i=1,nres + do j=1,3 + cref(j,i,1)=c(j,i) + cref(j,i+nres,1)=c(j,i+nres) + enddo + enddo + do i=1,2*nres + do j=1,3 + chomo(j,i,k)=c(j,i) + enddo + enddo + + return + end subroutine readpdb_template +!----------------------------------------------------------------------------- #endif !----------------------------------------------------------------------------- end module io_config