X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=inline;f=source%2Funres%2Fio_config.F90;h=5efb730665993a6840a01194bce6bc1ffae36ad4;hb=4db3a5f0b00eb4e0dc343054d4560d5e76f548f5;hp=2b1db84fa0c6217804eca5f5ee6135a1a37f2acf;hpb=b657faa3790a5b86aeeaa6f88aa849df89121b23;p=unres4.git diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90 index 2b1db84..5efb730 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 + sigmasctube,krad2,ract + integer :: ichir1,ichir2,ijunk,irdiff + character*3 string + character*80 temp1,mychar ! 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)) @@ -811,9 +813,9 @@ msc(:,:)=0.0d0 isc(:,:)=0.0d0 - allocate(msc(ntyp+1,5)) !(ntyp+1) - allocate(isc(ntyp+1,5)) !(ntyp+1) - allocate(restok(ntyp+1,5)) !(ntyp+1) + allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1) + allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1) + allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1) read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1) do i=1,ntyp_molec(1) @@ -826,6 +828,8 @@ dsc_inv(i)=1.0D0/dsc(i) endif enddo +! ip(1)=0.0001d0 +! isc(:,1)=0.0001d0 #endif read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2) do i=1,ntyp_molec(2) @@ -865,11 +869,18 @@ enddo enddo endif + + + + if (.not.allocated(ichargecat)) & + allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5))) + ichargecat(:)=0 + 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 +! ip(5)=0.2 ! isc(5)=0.2 read (iion,*) ncatprotparm allocate(catprm(ncatprotparm,4)) @@ -877,6 +888,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)) @@ -905,6 +933,36 @@ print *,liptranene(i) enddo close(iliptranpar) +! water parmaters entalphy + allocate(awaterenta(0:400)) + allocate(bwaterenta(0:400)) + allocate(cwaterenta(0:400)) + allocate(dwaterenta(0:400)) + allocate(awaterentro(0:400)) + allocate(bwaterentro(0:400)) + allocate(cwaterentro(0:400)) + allocate(dwaterentro(0:400)) + + read(iwaterwater,*) mychar + read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),& + cwaterenta(0),dwaterenta(0) + do i=1,398 + read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),& + cwaterenta(i),dwaterenta(i) + irdiff=int((ract-2.06d0)*50.0d0)+1 + if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff + enddo +! water parmaters entrophy + read(iwaterwater,*) mychar + read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),& + cwaterentro(0),dwaterentro(0) + do i=1,398 + read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),& + cwaterentro(i),dwaterentro(i) + irdiff=int((ract-2.06d0)*50.0d0)+1 + if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff + enddo + #ifdef CRYST_THETA ! @@ -1017,6 +1075,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 +1268,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 @@ -1436,6 +1518,41 @@ enddo endif enddo +#ifdef SC_END + allocate(nterm_scend(2,ntyp)) + allocate(arotam_end(0:6,2,ntyp)) + nterm_scend=0 + arotam_end=0.0d0 + read (irotam_end,*) ijunk +!c write (iout,*) "ijunk",ijunk + do i=1,ntyp + if (i.eq.10) cycle + do j=1,2 + read (irotam_end,'(a)') + read (irotam_end,*) nterm_scend(j,i) +!c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i) + do k=0,nterm_scend(j,i) + read (irotam_end,*) ijunk,arotam_end(k,j,i) +!c write (iout,*) "k",k," arotam",arotam_end(k,j,i) + enddo + enddo + enddo +!c lprint=.true. + if (lprint) then + write (iout,'(a)') & + "Parameters of the local potentials of sidechain ends" + do i=1,ntyp + write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')& + restyp(i,1),restyp(i,1) + do j=0,max0(nterm_scend(1,i),nterm_scend(2,i)) + write (iout,'(i5,2f20.10)') & + j,arotam_end(j,1,i),arotam_end(j,2,i) + enddo + enddo + endif +!c lprint=.false. +#endif + !---------reading nucleic acid parameters for rotamers------------------- allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp) do i=1,ntyp_molec(2) @@ -1501,6 +1618,455 @@ #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)) + bnew1=0.0d0 + bnew2=0.0d0 + ccnew=0.0d0 + ddnew=0.0d0 + e0new=0.0d0 + eenew=0.0d0 + bnew1tor=0.0d0 + bnew2tor=0.0d0 + ccnewtor=0.0d0 + ddnewtor=0.0d0 + e0newtor=0.0d0 + eenewtor=0.0d0 + 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 +2103,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 +2289,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 @@ -1889,6 +2569,7 @@ si=-1.0d0 do k=1,nterm_sccor(i,j) + print *,"test",i,j,k,l read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),& v2sccor(k,l,i,j) v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) @@ -1930,139 +2611,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 +2755,48 @@ 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 + read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp) + enddo + 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 @@ -2236,8 +2809,8 @@ 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) + dtail(1,i,j) = dtail(2,j,i) + dtail(2,i,j) = dtail(1,j,i) DO k = 1, 4 alphasur(k,i,j) = alphasur(k,j,i) wstate(k,i,j) = wstate(k,j,i) @@ -2258,10 +2831,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 +2929,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) @@ -2518,14 +3094,14 @@ allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2))) - + write (iout,*) "ESCBASEPARM" 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 + do j=1,ntyp_molec(2) ! 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) +! 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) @@ -2539,7 +3115,52 @@ read(isidep_scbase,*) & alphapol_scbase(i,j), & epsintab_scbase(i,j) + if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9 + if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9 + if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9 + if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9 + if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9 + if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9 + if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9 + if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9 + write(iout,*) & + 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) + write(iout,*) & + (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) + write(iout,*) & + dhead_scbasei(i,j), & + dhead_scbasej(i,j), & + rborn_scbasei(i,j),rborn_scbasej(i,j) + write(iout,*) & + (wdipdip_scbase(k,i,j),k=1,3), & + (wqdip_scbase(k,i,j),k=1,2) + write(iout,*) & + alphapol_scbase(i,j), & + epsintab_scbase(i,j) + END DO + j=4 + write(iout,*) & + 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) + write(iout,*) & + (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) + write(iout,*) & + dhead_scbasei(i,j), & + dhead_scbasej(i,j), & + rborn_scbasei(i,j),rborn_scbasej(i,j) + write(iout,*) & + (wdipdip_scbase(k,i,j),k=1,3), & + (wqdip_scbase(k,i,j),k=1,2) + write(iout,*) & + alphapol_scbase(i,j), & + epsintab_scbase(i,j) + END DO allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2))) @@ -2569,19 +3190,46 @@ allocate(chis_pepbase(ntyp_molec(2),2)) allocate(wdipdip_pepbase(3,ntyp_molec(2))) + write (iout,*) "EPEPBASEPARM" - do j=1,ntyp_molec(2)-1 ! without U then we will take T for U + do j=1,ntyp_molec(2) ! 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) + if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9 + if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9 + if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9 + if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9 + if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9 + if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9 + if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9 + if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9 + 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) + write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),& + chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2) + write(iout,*) & + (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), & + chis_pepbase(j,1),chis_pepbase(j,2) + write(iout,*) & + (wdipdip_pepbase(k,j),k=1,3) + END DO + j=4 + write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),& + chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2) + write(iout,*) & + (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), & + chis_pepbase(j,1),chis_pepbase(j,2) + write(iout,*) & + (wdipdip_pepbase(k,j),k=1,3) + allocate(aa_pepbase(ntyp_molec(2))) allocate(bb_pepbase(ntyp_molec(2))) @@ -2631,6 +3279,15 @@ read(isidep_scpho,*) & epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), & alphi_scpho(j) + if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9 + if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9 + if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9 + if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9 + if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9 + if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9 + if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9 + if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9 + END DO allocate(aa_scpho(ntyp_molec(1))) @@ -2762,19 +3419,304 @@ ! v1ss=0.0d0 ! v2ss=0.0d0 ! v3ss=0.0d0 - - if(me.eq.king) then - write (iout,'(/a)') "Disulfide bridge parameters:" - write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr - write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm - write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct - write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,& - ' v3ss:',v3ss +! MARTINI PARAMETER + allocate(ichargelipid(ntyp_molec(4))) + allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4))) + allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4))) + allocate(lip_bond(ntyp_molec(4),ntyp_molec(4))) + allocate(lip_eps(ntyp_molec(4),ntyp_molec(4))) + allocate(lip_sig(ntyp_molec(4),ntyp_molec(4))) + kjtokcal=0.2390057361 + krad=57.295779513 + !HERE THE MASS of MARTINI + write(*,*) "before MARTINI PARAM" + do i=1,ntyp_molec(4) + msc(i,4)=72.0d0 + mp(4)=0.0d0 + isc(i,4)=0.d0 + enddo + ip(4)=0.0 + msc(ntyp_molec(4)+1,4)=0.1d0 + !relative dielectric constant = 15 for implicit screening + k_coulomb_lip=332.0d0/15.0d0 + !kbond = 1250 kJ/(mol*nm*2) + kbondlip=1250.0d0*kjtokcal/100.0d0 + krad2=krad**2.0 + lip_angle_force=0.0d0 + if (DRY_MARTINI.gt.0) then + lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2 + lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2 + lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2 + lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2 + lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2 + lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2 + lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2 + lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2 + else + lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2 + lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2 + lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2 + lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2 + lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2 + lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2 + lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2 + lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2 endif - if (shield_mode.gt.0) then - pi=4.0D0*datan(1.0D0) -!C VSolvSphere the volume of solving sphere - print *,pi,"pi" + lip_angle_angle=0.0d0 + lip_angle_angle(3,12,12)=120.0/krad + lip_angle_angle(3,12,18)=180.0/krad + lip_angle_angle(3,18,16)=180.0/krad + lip_angle_angle(12,18,16)=180.0/krad + lip_angle_angle(18,16,18)=120.0/krad + lip_angle_angle(16,18,18)=180.0/krad + lip_angle_angle(12,18,18)=180.0/krad + lip_angle_angle(18,18,18)=180.0/krad + read(ilipbond,*) temp1 + do i=1,18 + read(ilipbond,*) temp1, lip_bond(i,1), & + lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), & + lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), & + lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), & + lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), & + lip_bond(i,18) + do j=1,18 + lip_bond(i,j)=lip_bond(i,j)*10 + enddo + enddo + + read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4)) + read(ilipnonbond,*) temp1 + do i=1,18 + read(ilipnonbond,*) temp1, lip_eps(i,1), & + lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), & + lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), & + lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), & + lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), & + lip_eps(i,18) +! write(*,*) i, lip_eps(i,18) + do j=1,18 + lip_eps(i,j)=lip_eps(i,j)*kjtokcal + enddo + enddo + read(ilipnonbond,*) temp1 + do i=1,18 + read(ilipnonbond,*) temp1,lip_sig(i,1), & + lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), & + lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), & + lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), & + lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), & + lip_sig(i,18) + do j=1,18 + lip_sig(i,j)=lip_sig(i,j)*10.0 + enddo + enddo + write(*,*) "after MARTINI PARAM" + +! Ions by Aga + + allocate(alphapolcat(ntyp,-1:ntyp_molec(5)),epsheadcat(ntyp,-1:ntyp_molec(5)),sig0headcat(ntyp,-1:ntyp_molec(5))) + allocate(alphapolcat2(ntyp,-1:ntyp_molec(5))) + allocate(sigiso1cat(ntyp,-1:ntyp_molec(5)),rborn1cat(ntyp,-1:ntyp_molec(5)),rborn2cat(ntyp,-1:ntyp_molec(5)),sigmap1cat(ntyp,-1:ntyp_molec(5))) + allocate(sigmap2cat(ntyp,-1:ntyp_molec(5)),sigiso2cat(ntyp,-1:ntyp_molec(5))) + allocate(chis1cat(ntyp,-1:ntyp_molec(5)),chis2cat(ntyp,-1:ntyp_molec(5)),wquadcat(ntyp,-1:ntyp_molec(5)),chipp1cat(ntyp,-1:ntyp_molec(5)),chipp2cat(ntyp,-1:ntyp_molec(5))) + allocate(epsintabcat(ntyp,-1:ntyp_molec(5))) + allocate(dtailcat(2,ntyp,-1:ntyp_molec(5))) + allocate(alphasurcat(4,ntyp,-1:ntyp_molec(5)),alphisocat(4,ntyp,-1:ntyp_molec(5))) + allocate(wqdipcat(2,ntyp,-1:ntyp_molec(5))) + allocate(wstatecat(4,ntyp,-1:ntyp_molec(5))) + allocate(dheadcat(2,2,ntyp,-1:ntyp_molec(5))) + allocate(nstatecat(ntyp,-1:ntyp_molec(5))) + allocate(debaykapcat(ntyp,-1:ntyp_molec(5))) + + if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,-1:ntyp1)) + if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,-1:ntyp1)) +! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp) + if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp) + if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp) + + + if (.not.allocated(ichargecat))& + allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5))) + write(*,*) "before ions",oldion + ichargecat(:)=0 + +! 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 + print *,ntyp_molec(5) + do i=-ntyp_molec(5),ntyp_molec(5) + read(iion,*) msc(i,5),restok(i,5),ichargecat(i) + print *,msc(i,5),restok(i,5) + enddo + ! ip(5)=0.2 + ! mp(5)=0.2 + pstok(5)=3.0 +!DIR$ NOUNROLL + do j=-1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals + if (j.eq.0) cycle + 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), & !6 + + (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&!12 +! chiscat(i,j),chiscat(j,i), & + chis1cat(i,j),chis2cat(i,j), & + + nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !19 !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),&!23 + dtailcat(1,i,j),dtailcat(2,i,j), & + epsheadcat(i,j),sig0headcat(i,j), &!27 +!wdipcat = w1 , w2 +! rborncat(i,j),rborncat(j,i),& + rborn1cat(i,j),rborn2cat(i,j),& + (wqdipcat(k,i,j),k=1,2), &!31 + alphapolcat(i,j),alphapolcat2(j,i), &!33 + (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j) + + if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(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,-1:ntyp_molec(5)),& + bb_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5))) + do i=1,ntyp + do j=-1,ntyp_molec(5) + if (j.eq.0) cycle + 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)-1 + 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(i,j) + write (iout,*) 'chip2= ', chipp2cat(i,j) + write (iout,*) 'alphasur1= ', alphasurcat(1,i,j) + write (iout,*) 'alphasur2= ', alphasurcat(2,i,j) + write (iout,*) 'alphasur3= ', alphasurcat(3,i,j) + write (iout,*) 'alphasur4= ', alphasurcat(4,i,j) + write (iout,*) 'sig1= ', sigmap1cat(i,j) + write (iout,*) 'sig2= ', sigmap2cat(i,j) + write (iout,*) 'chis1= ', chis1cat(i,j) + write (iout,*) 'chis1= ', chis2cat(i,j) + write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j) + write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j) + write (iout,*) 'dhead= ', dheadcat(1,1,i,j) + write (iout,*) 'dhead2= ', dheadcat(1,2,i,j) + write (iout,*) 'a1= ', rborn1cat(i,j) + write (iout,*) 'a2= ', rborn2cat(i,j) + write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i) + write (iout,*) 'alphapol1= ', alphapolcat(i,j) + write (iout,*) 'alphapol2= ', alphapolcat2(i,j) + write (iout,*) 'w1= ', wqdipcat(1,i,j) + write (iout,*) 'w2= ', wqdipcat(2,i,j) + write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,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 +! read number of Zn2+ +! here two denotes the Zn2+ and Cu2+ + write(iout,*) "before TRANPARM" + allocate(aomicattr(0:3,2)) + allocate(athetacattran(0:6,5,2)) + allocate(agamacattran(3,5,2)) + allocate(acatshiftdsc(5,2)) + allocate(bcatshiftdsc(5,2)) + allocate(demorsecat(5,2)) + allocate(alphamorsecat(5,2)) + allocate(x0catleft(5,2)) + allocate(x0catright(5,2)) + allocate(x0cattrans(5,2)) + allocate(ntrantyp(2)) + do i=1,1 ! currently only Zn2+ + + read(iiontran,*) ntrantyp(i) +!ntrantyp=4 +!| ao0 ao1 ao2 ao3 +!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5 +!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi +!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5 +!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5 + read(iiontran,*) (aomicattr(j,i),j=0,3) + do j=1,ntrantyp(i) + read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),& + (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),& + demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),& + x0cattrans(j,i) + enddo + enddo + if(me.eq.king) then + write (iout,'(/a)') "Disulfide bridge parameters:" + write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr + write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm + write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct + write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,& + ' v3ss:',v3ss + endif + if (shield_mode.gt.0) then + pi=4.0D0*datan(1.0D0) +!C VSolvSphere the volume of solving sphere + print *,pi,"pi" !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it !C there will be no distinction between proline peptide group and normal peptide !C group in case of shielding parameters @@ -2815,6 +3757,10 @@ 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" + go to 999 + 123 write(iout,*) "Error in transition metal parameters" 999 continue #ifdef MPI call MPI_Finalize(Ierror) @@ -2879,7 +3825,7 @@ character(len=3) :: seq,res,res2 character(len=5) :: atom character(len=80) :: card - real(kind=8),dimension(3,20) :: sccor + real(kind=8),dimension(3,40) :: sccor integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, integer :: isugar,molecprev,firstion character*1 :: sugar @@ -2952,6 +3898,7 @@ itype(ires_old,molecule)=ntyp1_molec(molecule) itype(ires_old-1,molecule)=ntyp1_molec(molecule) nres_molec(molecule)=nres_molec(molecule)+2 +! if (molecule.eq.4) ires=ires+2 ibeg=2 ! write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then @@ -2963,15 +3910,49 @@ ! iii=0 endif iii=0 + else if (card(:3).eq.'BRA') then + molecule=4 + ires=ires+1 + ires_old=ires+1 + itype(ires,molecule)=ntyp1_molec(molecule)-1 + nres_molec(molecule)=nres_molec(molecule)+1 + endif ! Read free energy 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 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then + if (card(14:16).eq.'LIP') then +! reading lipid + if (ibeg.eq.1) then + molecule=4 + ires=ires+1 + nres_molec(molecule)=nres_molec(molecule)+1 + itype(ires,molecule)=ntyp1_molec(molecule) + ibeg=0 + endif + if (ibeg.eq.2) then + ibeg=0 + ires=ires+2 + endif + + molecule=4 + nres_molec(molecule)=nres_molec(molecule)+1 + read (card(18:20),'(a3)') res + ires=ires+1 + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + if (UNRES_PDB) then! + itype(ires,molecule)=rescode(ires,res,0,molecule) + else + itype(ires,molecule)=rescode_lip(res,ires) + endif + else read (card(23:26),*) ires read (card(18:20),'(a3)') res ! write (iout,*) "ires",ires,ires-ishift+ishift1, @@ -2989,9 +3970,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 +3984,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 +4003,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 @@ -3046,11 +4027,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 @@ -3135,6 +4116,8 @@ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) endif endif + endif !LIP +! print *,"IONS",ions,card(1:6) else if ((ions).and.(card(1:6).eq.'HETATM')) then if (firstion.eq.0) then firstion=1 @@ -3144,24 +4127,26 @@ 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(1:2) read (card(18:20),'(a3)') res + if (atom(3:3).eq.'H') cycle if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.& (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') & - .or.(atom(1:2).eq.'K ')) & - then + .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.& + (atom(1:2).eq.'O ')) then ires=ires+1 + print *,"I have water" if (molecule.ne.5) molecprev=molecule molecule=5 nres_molec(molecule)=nres_molec(molecule)+1 print *,"HERE",nres_molec(molecule) - res=res(2:3)//' ' + if (res.ne.'WAT') 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 @@ -3235,7 +4220,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 @@ -3266,6 +4252,7 @@ ! print *,"molecule",molecule if ((itype(nres,1).ne.10)) then nres=nres+1 + nsup=nsup+1 if (molecule.eq.5) molecule=molecprev itype(nres,molecule)=ntyp1_molec(molecule) nres_molec(molecule)=nres_molec(molecule)+1 @@ -3282,8 +4269,8 @@ ! enddo ! else do j=1,3 - dcj=(c(j,nres-2)-c(j,nres-3))/2.0 - c(j,nres)=c(j,nres-1)+dcj + dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0 + c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj c(j,2*nres)=c(j,nres) enddo ! endif @@ -3350,7 +4337,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 @@ -3422,7 +4410,6 @@ do j=1,3 c_temporary(j,seqalingbegin)=c(j,i) c_temporary(j,seqalingbegin+nres)=c(j,i+nres) - enddo itype_temporary(seqalingbegin,k)=itype(i,k) print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin @@ -3435,6 +4422,9 @@ do i=1,2*nres do j=1,3 c(j,i)=c_temporary(j,i) + if (i.gt.nres) then + if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0 + endif enddo enddo do k=1,5 @@ -3443,6 +4433,24 @@ 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 j=1,3 + c(j,1)=c(j,2) + enddo + do i=2,nres-1 + itype(i,5)=itype(i+1,5) + do j=1,3 + c(j,i)=c(j,i+1) + enddo + 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 @@ -3478,7 +4486,7 @@ enddo endif -! print *,seqalingbegin,nres + print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -3514,13 +4522,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 @@ -3550,6 +4562,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) @@ -3559,16 +4578,16 @@ cou=1 write (iout,*) "symetr", symetr do i=1,nres - lll=lll+1 + lll=lll+1 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) - if (i.gt.1) then - if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then - chain_length=lll-1 - kkk=kkk+1 +! if (i.gt.1) then +! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then +! chain_length=lll-1 +! kkk=kkk+1 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) - lll=1 - endif - endif +! lll=1 +! endif +! endif do j=1,3 cref(j,i,cou)=c(j,i) cref(j,i+nres,cou)=c(j,i+nres) @@ -3578,13 +4597,6 @@ endif enddo enddo - write (iout,*) chain_length - if (chain_length.eq.0) chain_length=nres - do j=1,3 - chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1) - chain_rep(j,chain_length+nres,symetr) & - =chain_rep(j,chain_length+nres,1) - enddo ! diagnostic ! write (iout,*) "spraw lancuchy",chain_length,symetr ! do i=1,4 @@ -3599,36 +4611,36 @@ dc(j,0)=c(j,1) enddo - if (symetr.gt.1) then - call permut(symetr) - nperm=1 - do i=1,symetr - nperm=nperm*i - enddo - do i=1,nperm - write(iout,*) (tabperm(i,kkk),kkk=1,4) - enddo - do i=1,nperm - cou=0 - do kkk=1,symetr - icha=tabperm(i,kkk) - write (iout,*) i,icha - do lll=1,chain_length - cou=cou+1 - if (cou.le.nres) then - do j=1,3 - kupa=mod(lll,chain_length) - iprzes=(kkk-1)*chain_length+lll - if (kupa.eq.0) kupa=chain_length - write (iout,*) "kupa", kupa - cref(j,iprzes,i)=chain_rep(j,kupa,icha) - cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha) - enddo - endif - enddo - enddo - enddo - endif +! if (symetr.gt.1) then +! call permut(symetr) +! nperm=1 +! do i=1,symetr +! nperm=nperm*i +! enddo +! do i=1,nperm +! write(iout,*) (tabperm(i,kkk),kkk=1,4) +! enddo +! do i=1,nperm +! cou=0 +! do kkk=1,symetr +! icha=tabperm(i,kkk) +! write (iout,*) i,icha +! do lll=1,chain_length +! cou=cou+1 +! if (cou.le.nres) then +! do j=1,3 +! kupa=mod(lll,chain_length) +! iprzes=(kkk-1)*chain_length+lll +! if (kupa.eq.0) kupa=chain_length +! write (iout,*) "kupa", kupa +! cref(j,iprzes,i)=chain_rep(j,kupa,icha) +! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha) +! enddo +! endif +! enddo +! enddo +! enddo +! endif !-koniec robienia kopii ! diag do kkk=1,nperm @@ -3713,7 +4725,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 @@ -3735,6 +4748,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 @@ -3759,10 +4776,16 @@ 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 + fodson=index(controlcard,"FODSON").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')) @@ -3780,12 +4803,15 @@ ! 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) 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,"VNANO",velnanoconst,0.0d0) call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) @@ -3794,10 +4820,14 @@ 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 + call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0) + call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0) + call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0) + +! Lipidec parameters call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) if (lipthick.gt.0.0d0) then @@ -3935,6 +4965,7 @@ ! enddo buff_shield=1.0d0 endif + itime_mat=0 return end subroutine read_control !----------------------------------------------------------------------------- @@ -4033,7 +5064,7 @@ ! ! Read MD settings ! - use control_data, only: r_cut,rlamb,out1file + use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart use energy_data use geometry_data, only: pi use MPI_data @@ -4064,6 +5095,7 @@ call readi(controlcard,"NSTEP",n_timestep,1000000) call readi(controlcard,"NTWE",ntwe,100) call readi(controlcard,"NTWX",ntwx,1000) + call readi(controlcard,"NFOD",nfodstep,100) call reada(controlcard,"DT",d_time,1.0d-1) call reada(controlcard,"DVMAX",dvmax,2.0d1) call reada(controlcard,"DAMAX",damax,1.0d1) @@ -4079,6 +5111,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) @@ -4093,6 +5126,7 @@ print_compon = index(controlcard,"PRINT_COMPON").gt.0 rattle = index(controlcard,"RATTLE").gt.0 preminim=(index(controlcard,'PREMINIM').gt.0) + forceminim=(index(controlcard,'FORCEMINIM').gt.0) write (iout,*) "PREMINIM ",preminim dccart=(index(controlcard,'CART').gt.0) if (preminim) call read_minim @@ -4188,14 +5222,17 @@ enddo if(me.eq.king.or..not.out1file)then + do j=1,5 write (iout,'(/2a/)') & "Radii of site types and friction coefficients and std's of",& " stochastic forces of fully exposed sites" - write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1)) + write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j)) + do i=1,ntyp - write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& - gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1)) + write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),& + gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j)) enddo + enddo endif else if (tbf) then if(me.eq.king.or..not.out1file)then @@ -4603,6 +5640,10 @@ ! print *,"Processor",myrank," opened file ITHEP" call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',action='read') +#ifdef SC_END + call getenv_loc('ROTPAR_END',rotname_end) + open (irotam_end,file=rotname_end,status='old',action='read') +#endif ! print *,"Processor",myrank," opened file IROTAM" call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old',action='read') @@ -4639,7 +5680,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_TRAN',iontranname) + open (iiontran,file=iontranname,status='old',action='read') ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -4656,6 +5698,10 @@ open (ithep,file=thetname,status='old') call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old') +#ifdef SC_END + call getenv_loc('ROTPAR_END',rotname_end) + open (irotam_end,file=rotname_end,status='old') +#endif call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old') call getenv_loc('TORDPAR',tordname) @@ -4685,6 +5731,15 @@ 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') + call getenv_loc('IONPAR_TRAN',iontranname) + open (iiontran,file=iontranname,status='old') + call getenv_loc('WATWAT',iwaterwatername) + open (iwaterwater,file=iwaterwatername,status='old') + call getenv_loc('WATPROT',iwaterscname) + open (iwatersc,file=iwaterscname,status='old') + #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -4699,6 +5754,10 @@ open (ithep,file=thetname,status='old',action='read') call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',action='read') +#ifdef SC_END + call getenv_loc('ROTPAR_END',rotname_end) + open (irotam_end,file=rotname_end,status='old',action='read') +#endif call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old',action='read') call getenv_loc('TORDPAR',tordname) @@ -4740,10 +5799,22 @@ call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('LIPBOND',lipbondname) + open (ilipbond,file=lipbondname,status='old',action='read') + call getenv_loc('LIPNONBOND',lipnonbondname) + open (ilipnonbond,file=lipnonbondname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) 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') + call getenv_loc('IONPAR_TRAN',iontranname) + open (iiontran,file=iontranname,status='old',action='read') + call getenv_loc('WATWAT',iwaterwatername) + open (iwaterwater,file=iwaterwatername,status='old',action='read') + call getenv_loc('WATPROT',iwaterscname) + open (iwatersc,file=iwaterscname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) @@ -5347,5 +6418,339 @@ end subroutine write_stat_thread !----------------------------------------------------------------------------- #endif + 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,40) :: 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 +#if defined(WHAM_RUN) || defined(CLUSTER) + if (ires_old.lt.nres) then +#endif + itype(ires_old,1)=ntyp1 + iterter(ires_old)=1 +#if defined(WHAM_RUN) || defined(CLUSTER) + endif +#endif + 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 + write(2,*) "tutaj",ires,nres + 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