From b2ef3fdfb9fd19710ea49377832d6c3791f04c11 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 31 Mar 2021 08:40:41 +0200 Subject: [PATCH] critical bug fix for ions langvin and fix pdb output for wham and cluster --- source/cluster/io_clust.F90 | 7 ++++++- source/unres/MD.F90 | 8 ++++++-- source/unres/MREMD.F90 | 2 ++ source/unres/REMD.F90 | 8 ++++---- source/unres/energy.F90 | 8 ++++---- source/unres/geometry.F90 | 4 ++-- source/unres/io_config.F90 | 6 ++++-- source/wham/io_wham.F90 | 15 ++++++++++++--- 8 files changed, 40 insertions(+), 18 deletions(-) diff --git a/source/cluster/io_clust.F90 b/source/cluster/io_clust.F90 index ed4ccc3..3b3a134 100644 --- a/source/cluster/io_clust.F90 +++ b/source/cluster/io_clust.F90 @@ -1493,6 +1493,7 @@ call reada(weightcard,'WPEPBASE',wpepbase,1.0d0) call reada(weightcard,'WSCPHO',wscpho,0.0d0) call reada(weightcard,'WPEPPHO',wpeppho,0.0d0) + call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0) call reada(weightcard,"D0CM",d0cm,3.78d0) call reada(weightcard,"AKCM",akcm,15.1d0) @@ -1562,6 +1563,8 @@ weights(46)=wscbase weights(47)=wscpho weights(48)=wpeppho + weights(49)=wpeppho + weights(50)=wcatnucl write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,& @@ -1870,6 +1873,8 @@ write(iout,*)"po setup var" open (itube,file=tubename,status='old') call getenv('IONPAR',ionname) open (iion,file=ionname,status='old') + call getenv('IONPAR_NUCL',ionnuclname) + open (iionnucl,file=ionnuclname,status='old') #ifndef OLDSCP ! @@ -1946,7 +1951,7 @@ write(iout,*)"po setup var" write (ipdb,'(a)') 'TER' else ires=ires+1 - if (ires.eq.2) then + if ((ires.eq.2).and.(mnum.ne.5)) then do j=1,3 cbeg(j)=c(j,i-1) enddo diff --git a/source/unres/MD.F90 b/source/unres/MD.F90 index 2b412fd..3bb57e0 100644 --- a/source/unres/MD.F90 +++ b/source/unres/MD.F90 @@ -2058,6 +2058,7 @@ do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time +! write(iout,*) "adt",adt,"ads",adt2 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt d_t(j,inres)=d_t_old(j,inres)+adt @@ -5033,11 +5034,14 @@ !el ginvfric(2*nres,2*nres) !maxres2=2*maxres !el common /przechowalnia/ ginvfric - logical :: lprn = .false., checkmode = .false. + logical :: lprn, checkmode integer :: i,j,ind,k,nres2,nres6,mnum nres2=2*nres nres6=6*nres - + lprn=.false. + checkmode=.false. +! if (large) lprn=.true. +! if (large) checkmode=.true. if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres do i=0,nres2 diff --git a/source/unres/MREMD.F90 b/source/unres/MREMD.F90 index 8c4a4b5..003d67b 100644 --- a/source/unres/MREMD.F90 +++ b/source/unres/MREMD.F90 @@ -411,11 +411,13 @@ do i=nnt,nct-1 mnum=(molnum(i)) stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum)) +! write(iout,*) "stdforcp=",stdforcp(i),itype(i,mnum),i enddo do i=nnt,nct mnum=molnum(i) if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)& *dsqrt(gamsc(iabs(itype(i,mnum)),mnum)) +! write(iout,*) "stdforcsc=",stdforcsc(i),itype(i,mnum),i enddo endif call rescale_weights(t_bath) diff --git a/source/unres/REMD.F90 b/source/unres/REMD.F90 index fc0c9f1..699631c 100644 --- a/source/unres/REMD.F90 +++ b/source/unres/REMD.F90 @@ -627,11 +627,11 @@ ii = ind+m mnum=molnum(i) iti=itype(i,mnum) - if (mnum.eq.5) then - mscab=0.0 - else +! if (mnum.eq.5) then +! mscab=0.0 +! else mscab=msc(iabs(iti),mnum) - endif +! endif massvec(ii)=mscab if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index d235e01..ec3ffd4 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -11663,7 +11663,7 @@ endif endif endif -#define DEBUG +!#define DEBUG #ifdef DEBUG write (iout,*) "gradc gradx gloc" do i=1,nres @@ -11671,7 +11671,7 @@ i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) enddo #endif -#undef DEBUG +!#undef DEBUG #ifdef TIMING time_sumgradient=time_sumgradient+MPI_Wtime()-time01 #endif @@ -16730,13 +16730,13 @@ chip1=chip(itypi) #endif !#define DEBUG !el write (iout,*) "After sum_gradient" -!#ifdef DEBUG +#ifdef DEBUG write (iout,*) "After sum_gradient" do i=1,nres-1 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) enddo -!#endif +#endif !#undef DEBUG ! If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 index c89442f..fe8e89a 100644 --- a/source/unres/geometry.F90 +++ b/source/unres/geometry.F90 @@ -3398,7 +3398,7 @@ enddo endif ! Settind dE/ddnres-1 -#define DEBUG +!#define DEBUG #ifdef DEBUG j=1 write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), & @@ -3445,7 +3445,7 @@ write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), & gloc(ialph(i,1)+nside,icg),domega(j,3,i) #endif -#undef DEBUG +!#undef DEBUG enddo endif enddo diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90 index 73c18ab..0d92898 100644 --- a/source/unres/io_config.F90 +++ b/source/unres/io_config.F90 @@ -3293,6 +3293,8 @@ print *,msc(i,5),restok(i,5) enddo ip(5)=0.2 + ! mp(5)=0.2 + ! pstok(5)=1.0 !DIR$ NOUNROLL do j=1,ntyp_molec(5) do i=1,ntyp @@ -3935,8 +3937,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 diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90 index 68457dd..fee1e46 100644 --- a/source/wham/io_wham.F90 +++ b/source/wham/io_wham.F90 @@ -3908,8 +3908,9 @@ allocate(ww(max_eneW)) !-------------------------------------------------------------------------------- subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev) - use geometry_data, only:nres,c + use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum + use energy, only:boxshift ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -3924,7 +3925,7 @@ allocate(ww(max_eneW)) 'D','E','F','G','H','I','J','K','L','M','N','O',& 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid)) integer,dimension(nres) :: ica !(maxres) - real(kind=8) :: temp,efree,etot,entropy,rmsdev + real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj integer :: ii,i,j,iti,ires,iatom,ichain,mnum write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')& ii,temp,rmsdev @@ -3948,9 +3949,17 @@ allocate(ww(max_eneW)) ires=ires+1 iatom=iatom+1 ica(i)=iatom + if (mnum.ne.5) then write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,i),j=1,3) - if (iti.ne.10) then + else + xj=boxshift(c(1,i)-c(1,2),boxxsize) + yj=boxshift(c(2,i)-c(2,2),boxysize) + zj=boxshift(c(3,i)-c(3,2),boxzsize) + write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& + ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj + endif + if ((iti.ne.10).and.(mnum.ne.5)) then iatom=iatom+1 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,nres+i),j=1,3) -- 1.7.9.5