module ene_calc !----------------------------------------------------------------------------- use io_units use wham_data use control_data, only: tor_mode ! use geometry_data, only:nres,boxxsize,boxysize,boxzsize use energy_data ! use control_data, only:maxthetyp1 use energy, only:etotal,enerprint,rescale_weights #ifdef MPI use MPI_data ! include "mpif.h" ! include "COMMON.MPI" #endif implicit none !----------------------------------------------------------------------------- ! COMMON.ALLPARM ! common /allparm/ real(kind=8),dimension(:,:),allocatable :: ww_all !(max_ene,max_parm) ! max_eneW real(kind=8),dimension(:),allocatable :: vbldp0_all,akp_all !(max_parm) real(kind=8),dimension(:,:,:),allocatable :: vbldsc0_all,& aksc_all,abond0_all !(maxbondterm,ntyp,max_parm) real(kind=8),dimension(:,:),allocatable :: a0thet_all !(-ntyp:ntyp,max_parm) real(kind=8),dimension(:,:,:,:,:),allocatable :: athet_all,& bthet_all !(2,-ntyp:ntyp,-1:1,-1:1,max_parm) real(kind=8),dimension(:,:,:),allocatable :: polthet_all !(0:3,-ntyp:ntyp,max_parm) real(kind=8),dimension(:,:,:),allocatable :: gthet_all !(3,-ntyp:ntyp,max_parm) real(kind=8),dimension(:,:),allocatable :: theta0_all,& sig0_all,sigc0_all !(-ntyp:ntyp,max_parm) real(kind=8),dimension(:,:,:,:,:),allocatable :: aa0thet_all !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm) real(kind=8),dimension(:,:,:,:,:,:),allocatable :: aathet_all !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm) real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: bbthet_all,& ccthet_all,ddthet_all,eethet_all !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1, ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm) real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet_all1,& ggthet_all1,ffthet_all2,ggthet_all2 !(maxdouble,maxdouble,maxtheterm3, ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm) real(kind=8),dimension(:,:),allocatable :: dsc_all,dsc0_all !(ntyp1,max_parm) real(kind=8),dimension(:,:,:),allocatable :: bsc_all !(maxlob,ntyp,max_parm) real(kind=8),dimension(:,:,:,:),allocatable :: censc_all !(3,maxlob,-ntyp:ntyp,max_parm) real(kind=8),dimension(:,:,:,:,:),allocatable :: gaussc_all !(3,3,maxlob,-ntyp:ntyp,max_parm) real(kind=8),dimension(:,:,:),allocatable :: sc_parmin_all !(65,ntyp,max_parm) real(kind=8),dimension(:,:,:,:),allocatable :: v0_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm) real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_all,& v2_all !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) real(kind=8),dimension(:,:,:,:),allocatable :: vlor1_all,& vlor2_all,vlor3_all !(maxlor,maxtor,maxtor,max_parm) real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v1c_all,& v1s_all !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2c_all !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2s_all !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) real(kind=8),dimension(:,:,:),allocatable :: b1_all,b2_all !(2,-maxtor:maxtor,max_parm) real(kind=8),dimension(:,:,:,:),allocatable :: cc_all,dd_all,& ee_all !(2,2,-maxtor:maxtor,max_parm) real(kind=8),dimension(:,:,:,:),allocatable :: ctilde_all,& dtilde_all !(2,2,-maxtor:maxtor,max_parm) real(kind=8),dimension(:,:,:),allocatable :: b1tilde_all !(2,-maxtor:maxtor,max_parm) real(kind=8),dimension(:,:,:),allocatable :: app_all,bpp_all,& ael6_all,ael3_all !(2,2,max_parm) real(kind=8),dimension(:,:,:),allocatable :: aad_all,& bad_all !(ntyp,2,max_parm) real(kind=8),dimension(:,:,:),allocatable :: aa_all,bb_all,& augm_all,eps_all,sigma_all,r0_all,chi_all !(ntyp,ntyp,max_parm) real(kind=8),dimension(:,:),allocatable :: chip_all,alp_all !(ntyp,max_parm) real(kind=8),dimension(:),allocatable :: ebr_all,d0cm_all,& akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all !(max_parm) real(kind=8),dimension(:,:,:,:,:),allocatable :: v1sccor_all,& v2sccor_all !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm) integer,dimension(:,:),allocatable :: nlob_all !(ntyp1,max_parm) integer,dimension(:,:,:,:),allocatable :: nlor_all,& nterm_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm) integer,dimension(:,:,:,:,:),allocatable :: ntermd1_all,& ntermd2_all !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) integer,dimension(:,:),allocatable :: nbondterm_all !(ntyp,max_parm) integer,dimension(:,:),allocatable :: ithetyp_all !(-ntyp1:ntyp1,max_parm) integer,dimension(:),allocatable :: nthetyp_all,ntheterm_all,& ntheterm2_all,ntheterm3_all,nsingle_all,ndouble_all,& nntheterm_all !(max_parm) integer,dimension(:,:,:),allocatable :: nterm_sccor_all !(-ntyp:ntyp,-ntyp:ntyp,max_parm) !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- subroutine enecalc(islice,*) use names use control_data, only:indpdb use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,& vbld,rad2deg,dc_norm,dc,vbld_inv use io_base, only:gyrate!,briefout use geometry, only:int_from_cart1,returnbox use io_wham, only:pdboutW use io_database, only:opentmp use conform_compar, only:qwolynes,rmsnat ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" ! include "DIMENSIONS.FREE" #ifdef MPI ! use MPI_data include "mpif.h" ! include "COMMON.MPI" #endif ! include "COMMON.CHAIN" ! include "COMMON.IOUNITS" ! include "COMMON.PROTFILES" ! include "COMMON.NAMES" ! include "COMMON.VAR" ! include "COMMON.SBRIDGE" ! include "COMMON.GEO" ! include "COMMON.FFIELD" ! include "COMMON.ENEPS" ! include "COMMON.LOCAL" ! include "COMMON.WEIGHTS" ! include "COMMON.INTERACT" ! include "COMMON.FREE" ! include "COMMON.ENERGIES" ! include "COMMON.CONTROL" ! include "COMMON.TORCNSTR" ! implicit none #ifdef MPI integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) #endif character(len=64) :: nazwa character(len=80) :: bxname character(len=3) :: liczba !el real(kind=8) :: qwolynes !el external qwolynes integer :: errmsg_count,maxerrmsg_count=100 !el real(kind=8) :: rmsnat,gyrate !el external rmsnat,gyrate real(kind=8) :: tole=1.0d-1 integer i,itj,ii,iii,j,k,l,licz integer ir,ib,ipar,iparm integer iscor,islice real(kind=4) :: csingle(3,nres*2) real(kind=8) :: energ real(kind=8) :: temp !el integer ilen,iroof !el external ilen,iroof real(kind=8) :: energia(0:n_ene),rmsdev,efree,eini !el real(kind=8) :: energia(0:max_ene),rmsdev,efree,eini real(kind=8) :: fT(6),quot,quotl,kfacl,kfac=2.4d0,T0=3.0d2 real(kind=8) :: tt integer :: snk_p(MaxR,MaxT_h,nParmSet)!Max_parm) logical :: lerr character(len=64) :: bprotfile_temp ! integer :: rec integer,dimension(0:nprocs) :: scount_ !el real(kind=8) :: rmsnat rescale_mode=rescale_modeW call opentmp(islice,ientout,bprotfile_temp) iii=0 ii=0 !el ! iparm=1 errmsg_count=0 write (iout,*) "enecalc: nparmset ",nparmset #ifdef MPI do iparm=1,nParmSet do ib=1,nT_h(iparm) do i=1,nR(ib,iparm) snk_p(i,ib,iparm)=0 enddo enddo enddo do i=indstart(me1),indend(me1) write(iout,*)"enecalc_ i indstart",i,indstart(me1),indend(me1) #else do iparm=1,nParmSet do ib=1,nT_h(iparm) do i=1,nR(ib,iparm) snk(i,ib,iparm)=0 enddo enddo enddo do i=1,ntot write(iout,*)"enecalc_ i ntot",i,ntot #endif read(ientout,rec=i,err=101) & ((csingle(l,k),l=1,3),k=1,nres),& ((csingle(l,k+nres),l=1,3),k=nnt,nct),& nss,(ihpb(k),jhpb(k),k=1,nss),& eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar !el debug !write(iout,*)"co wczytuje" ! write(iout,*)((csingle(l,k),l=1,3),k=1,nres),& ! ((csingle(l,k+nres),l=1,3),k=nnt,nct),& ! nss,(ihpb(k),jhpb(k),k=1,nss),& ! eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar !el -------- !write(iout,*)"ipar",ib,ipar,1.0d0/(beta_h(ib,ipar)*1.987D-3) if (indpdb.gt.0) then do k=1,nres do l=1,3 c(l,k)=csingle(l,k) enddo enddo do k=nnt,nct do l=1,3 c(l,k+nres)=csingle(l,k+nres) enddo enddo call returnbox do k=1,nres do l=1,3 csingle(l,k)=c(l,k) enddo enddo do k=nnt,nct do l=1,3 csingle(l,k+nres)=c(l,k+nres) enddo enddo anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3) q(nQ+1,iii+1)=rmsnat(iii+1) endif q(nQ+2,iii+1)=gyrate(iii+1) ! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el ! fT=T0*beta_h(ib,ipar)*1.987D-3 ! ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)) ! EL start old rescale ! if (rescale_modeW.eq.1) then ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3) !#if defined(FUNCTH) ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3) ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0 !#elif defined(FUNCT) ! ft(6)=quot !#else ! ft(6)=1.0d0 !#endif ! quotl=1.0d0 ! kfacl=1.0d0 ! do l=1,5 ! quotl=quotl*quot ! kfacl=kfacl*kfac ! fT(l)=kfacl/(kfacl-1.0d0+quotl) ! enddo ! else if (rescale_modeW.eq.2) then ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3) !#if defined(FUNCTH) ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3) ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0 !#elif defined(FUNCT) ! ft(6)=quot !#else ! ft(6)=1.0d0 !#endif ! quotl=1.0d0 ! do l=1,5 ! quotl=quotl*quot ! fT(l)=1.12692801104297249644d0/ & ! dlog(dexp(quotl)+dexp(-quotl)) ! enddo ! else if (rescale_modeW.eq.0) then ! do l=1,5 ! fT(l)=1.0d0 ! enddo ! else ! write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",& ! rescale_modeW ! call flush(iout) ! return 1 ! endif !EL end old rescele ! write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0, ! & " kfac",kfac,"quot",quot," fT",fT #ifdef DEBUG write(iout,*)"weights" write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,& wtor_d,wsccor,wbond #endif do j=1,2*nres do k=1,3 c(k,j)=csingle(k,j) enddo enddo call int_from_cart1(.false.) ii=ii+1 ! call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3)) do iparm=1,nparmset #ifdef DEBUG write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3) write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,& wtor_d,wsccor,wbond #endif call restore_parm(iparm) call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3)) #ifdef DEBUG write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3) write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,& wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,& wtor_d,wsccor,wbond #endif ! call etotal(energia(0),fT) call etotal(energia(0)) !write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar)) !do k=1,2*nres+2 !write(iout,*)k,"c=",(c(l,k),l=1,3) !write(iout,*)k,"dc=",(dc(l,k),l=1,3) !write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3) !enddo !do k=1,nres*2 !write(iout,*)k,"vbld=",vbld(k) !write(iout,*)k,"vbld_inv=",vbld_inv(k) !enddo !write(iout,*)"energia",(energia(j),j=0,n_ene) !write(iout,*)"enerprint tuz po call etotal" call enerprint(energia(0)) #ifdef DEBUG write (iout,*) "Conformation",i write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) ! call enerprint(energia(0),fT) call enerprint(energia(0)) write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49) write (iout,*) "ftors",ftors !el call briefout(i,energia(0)) temp=1.0d0/(beta_h(ib,ipar)*1.987D-3) write (iout,*) "temp", temp call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0) #endif if (energia(0).ge.1.0d20) then write (iout,*) "NaNs detected in some of the energy",& " components for conformation",ii+1 write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" ! call intout ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) "The components of the energy are:" ! call enerprint(energia(0),fT) call enerprint(energia(0)) write (iout,*) & "This conformation WILL NOT be added to the database." call flush(iout) goto 121 else #ifdef DEBUG if (ipar.eq.iparm) write (iout,*) i,iparm,& 1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0) #endif if (ipar.eq.iparm .and. einicheck.gt.0 .and. & dabs(eini-energia(0)).gt.tole) then if (errmsg_count.le.maxerrmsg_count) then write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') & "Warning: energy differs remarkably from ",& " the value read in: ",energia(0),eini," point",& iii+1,indstart(me1)+iii," T",& 1.0d0/(1.987D-3*beta_h(ib,ipar)) ! call intout call pdboutW(indstart(me1)+iii,& 1.0d0/(1.987D-3*beta_h(ib,ipar)),& energia(0),eini,0.0d0,0.0d0) write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" ! call intout ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) call enerprint(energia(0)) ! call enerprint(energia(0),fT) errmsg_count=errmsg_count+1 if (errmsg_count.gt.maxerrmsg_count) & write (iout,*) "Too many warning messages" if (einicheck.gt.1) then write (iout,*) "Calculation stopped." call flush(iout) #ifdef MPI call MPI_Abort(WHAM_COMM,IERROR,ERRCODE) #endif call flush(iout) return 1 endif endif endif potE(iii+1,iparm)=energia(0) do k=1,49 enetb(k,iii+1,iparm)=energia(k) enddo ! write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21) ! write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21) #ifdef DEBUG write (iout,'(2i5,f10.1,3e15.5)') i,iii,& 1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree ! call enerprint(energia(0),fT) #endif #ifdef DEBUG write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss) write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ) write (iout,'(f10.5,i10)') rmsdev,iscor ! call enerprint(energia(0),fT) call enerprint(energia(0)) call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) #endif endif enddo ! iparm iii=iii+1 if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0) write (ientout,rec=iii) & ((csingle(l,k),l=1,3),k=1,nres),& ((csingle(l,k+nres),l=1,3),k=nnt,nct),& nss,(ihpb(k),jhpb(k),k=1,nss),& potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar ! write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree #ifdef MPI if (separate_parset) then snk_p(iR,ib,1)=snk_p(iR,ib,1)+1 else snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1 endif ! write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar, ! & " snk",snk_p(iR,ib,ipar) #else snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1 #endif 121 continue enddo #ifdef MPI scount(me)=iii write (iout,*) "Me",me," scount",scount(me) call flush(iout) ! Master gathers updated numbers of conformations written by all procs. call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, & MPI_INTEGER, WHAM_COMM, IERROR) indstart(0)=1 indend(0)=scount_(0) do i=1, Nprocs-1 indstart(i)=indend(i-1)+1 indend(i)=indstart(i)+scount_(i)-1 enddo write (iout,*) write (iout,*) "Revised conformation counts" do i=0,nprocs1-1 write (iout,'(a,i5,a,i7,a,i7,a,i7)') & "Processor",i," indstart",indstart(i),& " indend",indend(i)," count",scount_(i) enddo call flush(iout) call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),& MaxR*MaxT_h*nParmSet,& MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR) #endif stot(islice)=0 do iparm=1,nParmSet do ib=1,nT_h(iparm) do i=1,nR(ib,iparm) stot(islice)=stot(islice)+snk(i,ib,iparm,islice) enddo enddo enddo write (iout,*) "Revised SNK" do iparm=1,nParmSet do ib=1,nT_h(iparm) write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') & iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),& (snk(i,ib,iparm,islice),i=1,nR(ib,iparm)) write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm)) enddo enddo write (iout,'("Total",i10)') stot(islice) call flush(iout) do i=0,nprocs scount(i)=scount_(i) enddo return 101 write (iout,*) "Error in scratchfile." call flush(iout) !el#undef DEBUG return 1 end subroutine enecalc !------------------------------------------------------------------------------ logical function conf_check(ii,iprint) use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld use geometry, only:int_from_cart1 ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" ! include "DIMENSIONS.FREE" #ifdef MPI ! use MPI_data include "mpif.h" ! include "COMMON.MPI" #endif ! include "COMMON.CHAIN" ! include "COMMON.IOUNITS" ! include "COMMON.PROTFILES" ! include "COMMON.NAMES" ! include "COMMON.VAR" ! include "COMMON.SBRIDGE" ! include "COMMON.GEO" ! include "COMMON.FFIELD" ! include "COMMON.ENEPS" ! include "COMMON.LOCAL" ! include "COMMON.WEIGHTS" ! include "COMMON.INTERACT" ! include "COMMON.FREE" ! include "COMMON.ENERGIES" ! include "COMMON.CONTROL" ! include "COMMON.TORCNSTR" ! implicit none #ifdef MPI integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) #endif integer :: j,k,l,ii,itj,iprint,mnum,mnum1 if (.not.check_conf) then conf_check=.true. return endif call int_from_cart1(.false.) do j=nnt+1,nct mnum=molnum(j) mnum1=molnum(j-1) if (mnum.ne.1) cycle if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. & (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then if (iprint.gt.0) & write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),& " for conformation",ii if (iprint.gt.1) then write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) endif if (iprint.gt.0) write (iout,*) & "This conformation WILL NOT be added to the database." conf_check=.false. return endif enddo do j=nnt,nct mnum=molnum(j) if (mnum.gt.3) cycle itj=itype(j,mnum) if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. & (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then if (iprint.gt.0) & write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),& " for conformation",ii if (iprint.gt.1) then write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) endif if (iprint.gt.0) write (iout,*) & "This conformation WILL NOT be added to the database." conf_check=.false. return endif enddo do j=3,nres if (theta(j).le.0.0d0) then if (iprint.gt.0) & write (iout,*) "Zero theta angle(s) in conformation",ii if (iprint.gt.1) then write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) endif if (iprint.gt.0) write (iout,*) & "This conformation WILL NOT be added to the database." conf_check=.false. return endif if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad enddo conf_check=.true. ! write (iout,*) "conf_check passed",ii return end function conf_check !----------------------------------------------------------------------------- ! store_parm.F !----------------------------------------------------------------------------- subroutine store_parm(iparm) ! ! Store parameters of set IPARM ! valence angles and the side chains and energy parameters. ! ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.FREE' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.TORSION' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' ! include 'COMMON.SBRIDGE' ! include 'COMMON.SCROT' ! include 'COMMON.SCCOR' ! include 'COMMON.ALLPARM' integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii call alloc_enecalc_arrays(iparm) !el allocate(ww_all(n_ene,iparm)) ! Store weights ww_all(1,iparm)=wsc ww_all(2,iparm)=wscp ww_all(3,iparm)=welec ww_all(4,iparm)=wcorr ww_all(5,iparm)=wcorr5 ww_all(6,iparm)=wcorr6 ww_all(7,iparm)=wel_loc ww_all(8,iparm)=wturn3 ww_all(9,iparm)=wturn4 ww_all(10,iparm)=wturn6 ww_all(11,iparm)=wang ww_all(12,iparm)=wscloc ww_all(13,iparm)=wtor ww_all(14,iparm)=wtor_d ww_all(15,iparm)=wstrain ww_all(16,iparm)=wvdwpp ww_all(17,iparm)=wbond ww_all(19,iparm)=wsccor ww_all(42,iparm)=wcatprot ww_all(41,iparm)=wcatcat ! Store bond parameters vbldp0_all(iparm)=vbldp0 akp_all(iparm)=akp do i=1,ntyp nbondterm_all(i,iparm)=nbondterm(i) do j=1,nbondterm(i) vbldsc0_all(j,i,iparm)=vbldsc0(j,i) aksc_all(j,i,iparm)=aksc(j,i) abond0_all(j,i,iparm)=abond0(j,i) enddo enddo ! Store bond angle parameters #ifdef CRYST_THETA do i=-ntyp,ntyp a0thet_all(i,iparm)=a0thet(i) do ichir1=-1,1 do ichir2=-1,1 do j=1,2 athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2) bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2) enddo enddo enddo do j=0,3 polthet_all(j,i,iparm)=polthet(j,i) enddo do j=1,3 gthet_all(j,i,iparm)=gthet(j,i) enddo theta0_all(i,iparm)=theta0(i) sig0_all(i,iparm)=sig0(i) sigc0_all(i,iparm)=sigc0(i) enddo #else IF (tor_mode.eq.0) THEN nthetyp_all(iparm)=nthetyp ntheterm_all(iparm)=ntheterm ntheterm2_all(iparm)=ntheterm2 ntheterm3_all(iparm)=ntheterm3 nsingle_all(iparm)=nsingle ndouble_all(iparm)=ndouble nntheterm_all(iparm)=nntheterm do i=-ntyp,ntyp ithetyp_all(i,iparm)=ithetyp(i) enddo do iblock=1,2 do i=-nthetyp-1,nthetyp+1 do j=-nthetyp-1,nthetyp+1 do k=-nthetyp-1,nthetyp+1 aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock) do l=1,ntheterm aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock) enddo do l=1,ntheterm2 do m=1,nsingle bbthet_all(m,l,i,j,k,iblock,iparm)= & bbthet(m,l,i,j,k,iblock) ccthet_all(m,l,i,j,k,iblock,iparm)= & ccthet(m,l,i,j,k,iblock) ddthet_all(m,l,i,j,k,iblock,iparm)= & ddthet(m,l,i,j,k,iblock) eethet_all(m,l,i,j,k,iblock,iparm)= & eethet(m,l,i,j,k,iblock) enddo enddo do l=1,ntheterm3 do m=1,ndouble do mm=1,ndouble if (iblock.eq.1) then ffthet_all1(mm,m,l,i,j,k,iparm)=& ffthet(mm,m,l,i,j,k,iblock) ggthet_all1(mm,m,l,i,j,k,iparm)=& ggthet(mm,m,l,i,j,k,iblock) else ffthet_all2(mm,m,l,i,j,k,iparm)=& ffthet(mm,m,l,i,j,k,iblock) ggthet_all2(mm,m,l,i,j,k,iparm)=& ggthet(mm,m,l,i,j,k,iblock) endif enddo enddo enddo enddo enddo enddo enddo ELSE write(iout,*) "Need storing parameters" ENDIF #endif #ifdef CRYST_SC ! Store the sidechain rotamer parameters do i=-ntyp,ntyp iii=iabs(i) !! write (iout,*) i,"storeparm1" if (i.eq.0) cycle nlob_all(iii,iparm)=nlob(iii) do j=1,nlob(iii) bsc_all(j,iii,iparm)=bsc(j,iii) do k=1,3 censc_all(k,j,i,iparm)=censc(k,j,i) enddo do k=1,3 do l=1,3 gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i) enddo enddo enddo enddo #else do i=1,ntyp do j=1,65 sc_parmin_all(j,i,iparm)=sc_parmin(j,i) enddo enddo #endif IF (tor_mode.eq.0) THEN ! Store the torsional parameters do iblock=1,2 do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 v0_all(i,j,iblock,iparm)=v0(i,j,iblock) nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock) nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock) do k=1,nterm(i,j,iblock) v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock) v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock) enddo do k=1,nlor(i,j,iblock) vlor1_all(k,i,j,iparm)=vlor1(k,i,j) vlor2_all(k,i,j,iparm)=vlor2(k,i,j) vlor3_all(k,i,j,iparm)=vlor3(k,i,j) enddo enddo enddo enddo ! Store the double torsional parameters do iblock=1,2 do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 do k=-ntortyp+1,ntortyp-1 ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock) ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock) do l=1,ntermd_1(i,j,k,iblock) v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock) v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock) v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock) v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock) enddo do l=1,ntermd_2(i,j,k,iblock) do m=1,ntermd_2(i,j,k,iblock) v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock) enddo enddo enddo enddo enddo enddo ! Store parameters of the cumulants do i=-nloctyp,nloctyp do j=1,2 b1_all(j,i,iparm)=b1(j,i) b1tilde_all(j,i,iparm)=b1tilde(j,i) b2_all(j,i,iparm)=b2(j,i) enddo do j=1,2 do k=1,2 cc_all(k,j,i,iparm)=cc(k,j,i) ctilde_all(k,j,i,iparm)=ctilde(k,j,i) dd_all(k,j,i,iparm)=dd(k,j,i) dtilde_all(k,j,i,iparm)=dtilde(k,j,i) ee_all(k,j,i,iparm)=ee(k,j,i) enddo enddo enddo ELSE write(iout,*) "NEED storing parameters" ENDIF ! Store the parameters of electrostatic interactions do i=1,2 do j=1,2 app_all(j,i,iparm)=app(j,i) bpp_all(j,i,iparm)=bpp(j,i) ael6_all(j,i,iparm)=ael6(j,i) ael3_all(j,i,iparm)=ael3(j,i) enddo enddo ! Store sidechain parameters do i=1,ntyp do j=1,ntyp aa_all(j,i,iparm)=aa_aq(j,i) bb_all(j,i,iparm)=bb_aq(j,i) r0_all(j,i,iparm)=r0(j,i) sigma_all(j,i,iparm)=sigma(j,i) chi_all(j,i,iparm)=chi(j,i) augm_all(j,i,iparm)=augm(j,i) eps_all(j,i,iparm)=eps(j,i) enddo enddo do i=1,ntyp chip_all(i,iparm)=chip(i) alp_all(i,iparm)=alp(i) enddo ! Store the SCp parameters do i=1,ntyp do j=1,2 aad_all(i,j,iparm)=aad(i,j) bad_all(i,j,iparm)=bad(i,j) enddo enddo ! Store disulfide-bond parameters ebr_all(iparm)=ebr d0cm_all(iparm)=d0cm akcm_all(iparm)=akcm akth_all(iparm)=akth akct_all(iparm)=akct v1ss_all(iparm)=v1ss v2ss_all(iparm)=v2ss v3ss_all(iparm)=v3ss ! Store SC-backbone correlation parameters do i=-nsccortyp,nsccortyp do j=-nsccortyp,nsccortyp nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i) ! do i=1,20 ! do j=1,20 do l=1,3 do k=1,nterm_sccor(j,i) v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i) v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i) enddo enddo enddo enddo write(iout,*)"end of store_parm" return end subroutine store_parm !-------------------------------------------------------------------------- subroutine restore_parm(iparm) ! ! Store parameters of set IPARM ! valence angles and the side chains and energy parameters. ! ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.FREE' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.TORSION' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' ! include 'COMMON.SBRIDGE' ! include 'COMMON.SCROT' ! include 'COMMON.SCCOR' ! include 'COMMON.ALLPARM' integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii ! Restore weights wsc=ww_all(1,iparm) wscp=ww_all(2,iparm) welec=ww_all(3,iparm) wcorr=ww_all(4,iparm) wcorr5=ww_all(5,iparm) wcorr6=ww_all(6,iparm) wel_loc=ww_all(7,iparm) wturn3=ww_all(8,iparm) wturn4=ww_all(9,iparm) wturn6=ww_all(10,iparm) wang=ww_all(11,iparm) wscloc=ww_all(12,iparm) wtor=ww_all(13,iparm) wtor_d=ww_all(14,iparm) wstrain=ww_all(15,iparm) wvdwpp=ww_all(16,iparm) wbond=ww_all(17,iparm) wsccor=ww_all(19,iparm) wcatprot=ww_all(42,iparm) wcatcat=ww_all(41,iparm) ! Restore bond parameters vbldp0=vbldp0_all(iparm) akp=akp_all(iparm) do i=1,ntyp nbondterm(i)=nbondterm_all(i,iparm) do j=1,nbondterm(i) vbldsc0(j,i)=vbldsc0_all(j,i,iparm) aksc(j,i)=aksc_all(j,i,iparm) abond0(j,i)=abond0_all(j,i,iparm) enddo enddo ! Restore bond angle parameters #ifdef CRYST_THETA do i=-ntyp,ntyp a0thet(i)=a0thet_all(i,iparm) do ichir1=-1,1 do ichir2=-1,1 do j=1,2 athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm) bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm) enddo enddo enddo do j=0,3 polthet(j,i)=polthet_all(j,i,iparm) enddo do j=1,3 gthet(j,i)=gthet_all(j,i,iparm) enddo theta0(i)=theta0_all(i,iparm) sig0(i)=sig0_all(i,iparm) sigc0(i)=sigc0_all(i,iparm) enddo #else IF (tor_mode.eq.0) THEN nthetyp=nthetyp_all(iparm) ntheterm=ntheterm_all(iparm) ntheterm2=ntheterm2_all(iparm) ntheterm3=ntheterm3_all(iparm) nsingle=nsingle_all(iparm) ndouble=ndouble_all(iparm) nntheterm=nntheterm_all(iparm) do i=-ntyp,ntyp ithetyp(i)=ithetyp_all(i,iparm) enddo do iblock=1,2 do i=-nthetyp-1,nthetyp+1 do j=-nthetyp-1,nthetyp+1 do k=-nthetyp-1,nthetyp+1 aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm) do l=1,ntheterm aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm) enddo do l=1,ntheterm2 do m=1,nsingle bbthet(m,l,i,j,k,iblock)= & bbthet_all(m,l,i,j,k,iblock,iparm) ccthet(m,l,i,j,k,iblock)= & ccthet_all(m,l,i,j,k,iblock,iparm) ddthet(m,l,i,j,k,iblock)= & ddthet_all(m,l,i,j,k,iblock,iparm) eethet(m,l,i,j,k,iblock)= & eethet_all(m,l,i,j,k,iblock,iparm) enddo enddo do l=1,ntheterm3 do m=1,ndouble do mm=1,ndouble if (iblock.eq.1) then ffthet(mm,m,l,i,j,k,iblock)= & ffthet_all1(mm,m,l,i,j,k,iparm) ggthet(mm,m,l,i,j,k,iblock)= & ggthet_all1(mm,m,l,i,j,k,iparm) else ffthet(mm,m,l,i,j,k,iblock)= & ffthet_all2(mm,m,l,i,j,k,iparm) ggthet(mm,m,l,i,j,k,iblock)= & ggthet_all2(mm,m,l,i,j,k,iparm) endif enddo enddo enddo enddo enddo enddo enddo ELSE write (iout,*) "Need storing parameters" ENDIF #endif ! Restore the sidechain rotamer parameters #ifdef CRYST_SC do i=-ntyp,ntyp if (i.eq.0) cycle iii=iabs(i) nlob(iii)=nlob_all(iii,iparm) do j=1,nlob(iii) bsc(j,iii)=bsc_all(j,iii,iparm) do k=1,3 censc(k,j,i)=censc_all(k,j,i,iparm) enddo do k=1,3 do l=1,3 gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm) enddo enddo enddo enddo #else do i=1,ntyp do j=1,65 sc_parmin(j,i)=sc_parmin_all(j,i,iparm) enddo enddo #endif IF (tor_mode.eq.0) THEN ! Restore the torsional parameters do iblock=1,2 do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 v0(i,j,iblock)=v0_all(i,j,iblock,iparm) nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm) nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm) do k=1,nterm(i,j,iblock) v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm) v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm) enddo do k=1,nlor(i,j,iblock) vlor1(k,i,j)=vlor1_all(k,i,j,iparm) vlor2(k,i,j)=vlor2_all(k,i,j,iparm) vlor3(k,i,j)=vlor3_all(k,i,j,iparm) enddo enddo enddo enddo ! Restore the double torsional parameters do iblock=1,2 do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 do k=-ntortyp+1,ntortyp-1 ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm) ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm) do l=1,ntermd_1(i,j,k,iblock) v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm) v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm) v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm) v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm) enddo do l=1,ntermd_2(i,j,k,iblock) do m=1,ntermd_2(i,j,k,iblock) v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm) enddo enddo enddo enddo enddo enddo ! Restore parameters of the cumulants do i=-nloctyp,nloctyp do j=1,2 b1(j,i)=b1_all(j,i,iparm) b1tilde(j,i)=b1tilde_all(j,i,iparm) b2(j,i)=b2_all(j,i,iparm) enddo do j=1,2 do k=1,2 cc(k,j,i)=cc_all(k,j,i,iparm) ctilde(k,j,i)=ctilde_all(k,j,i,iparm) dd(k,j,i)=dd_all(k,j,i,iparm) dtilde(k,j,i)=dtilde_all(k,j,i,iparm) ee(k,j,i)=ee_all(k,j,i,iparm) enddo enddo enddo ELSE write(iout,*) "need storing parameters" ENDIF ! Restore the parameters of electrostatic interactions do i=1,2 do j=1,2 app(j,i)=app_all(j,i,iparm) bpp(j,i)=bpp_all(j,i,iparm) ael6(j,i)=ael6_all(j,i,iparm) ael3(j,i)=ael3_all(j,i,iparm) enddo enddo ! Restore sidechain parameters do i=1,ntyp do j=1,ntyp aa_aq(j,i)=aa_all(j,i,iparm) bb_aq(j,i)=bb_all(j,i,iparm) r0(j,i)=r0_all(j,i,iparm) sigma(j,i)=sigma_all(j,i,iparm) chi(j,i)=chi_all(j,i,iparm) augm(j,i)=augm_all(j,i,iparm) eps(j,i)=eps_all(j,i,iparm) enddo enddo do i=1,ntyp chip(i)=chip_all(i,iparm) alp(i)=alp_all(i,iparm) enddo ! Restore the SCp parameters do i=1,ntyp do j=1,2 aad(i,j)=aad_all(i,j,iparm) bad(i,j)=bad_all(i,j,iparm) enddo enddo ! Restore disulfide-bond parameters ebr=ebr_all(iparm) d0cm=d0cm_all(iparm) akcm=akcm_all(iparm) akth=akth_all(iparm) akct=akct_all(iparm) v1ss=v1ss_all(iparm) v2ss=v2ss_all(iparm) v3ss=v3ss_all(iparm) ! Restore SC-backbone correlation parameters do i=-nsccortyp,nsccortyp do j=-nsccortyp,nsccortyp nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm) do l=1,3 do k=1,nterm_sccor(j,i) v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm) v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm) enddo enddo enddo enddo return end subroutine restore_parm !-------------------------------------------------------------------------- ! make_ensemble1.F !-------------------------------------------------------------------------- subroutine make_ensembles(islice,*) ! construct the conformational ensembles at REMD temperatures use geometry_data, only:c use io_base, only:ilen use io_wham, only:pdboutW use geometry, only:returnbox ! implicit none ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" ! include "DIMENSIONS.FREE" #ifdef MPI include "mpif.h" ! include "COMMON.MPI" integer :: ierror,errcode,status(MPI_STATUS_SIZE) #endif ! include "COMMON.IOUNITS" ! include "COMMON.CONTROL" ! include "COMMON.FREE" ! include "COMMON.ENERGIES" ! include "COMMON.FFIELD" ! include "COMMON.INTERACT" ! include "COMMON.SBRIDGE" ! include "COMMON.CHAIN" ! include "COMMON.PROTFILES" ! include "COMMON.PROT" real(kind=4) :: csingle(3,nres*2) real(kind=8),dimension(6) :: fT,fTprim,fTbis real(kind=8) :: quot,quotl1,quotl,kfacl,& eprim,ebis,temper,kfac=2.4d0,T0=300.0d0 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,& eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, & ecation_prot, ecationcation integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist real(kind=8) :: qfree,sumprob,eini,efree,rmsdev character(len=80) :: bxname character(len=2) :: licz1,licz2 character(len=3) :: licz3,licz4 character(len=5) :: ctemper !el integer ilen !el external ilen real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr) real(kind=8) :: enepot(MaxStr) integer :: iperm(MaxStr) integer :: islice integer,dimension(0:nprocs) :: scount_ #ifdef MPI if (me.eq.Master) then #endif write (licz2,'(bz,i2.2)') islice if (nslice.eq.1) then if (.not.separate_parset) then bxname = prefix(:ilen(prefix))//".bx" else write (licz3,'(bz,i3.3)') myparm bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx" endif else if (.not.separate_parset) then bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx" else write (licz3,'(bz,i3.3)') myparm bxname = prefix(:ilen(prefix))//"par_"//licz3// & "_slice_"//licz2//".bx" endif endif open (ientout,file=bxname,status="unknown",& form="unformatted",access="direct",recl=lenrec1) #ifdef MPI endif #endif do iparm=1,iparm if (iparm.ne.iparmprint) exit call restore_parm(iparm) do ib=1,nT_h(iparm) #ifdef DEBUG write (iout,*) "iparm",iparm," ib",ib #endif temper=1.0d0/(beta_h(ib,iparm)*1.987D-3) ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) ! quotl=1.0d0 ! kfacl=1.0d0 ! do l=1,5 ! quotl1=quotl ! quotl=quotl*quot ! kfacl=kfacl*kfac ! fT(l)=kfacl/(kfacl-1.0d0+quotl) ! enddo !el old rescale weights ! ! if (rescale_mode.eq.1) then ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) #if defined(FUNCTH) tt=1.0d0/(beta_h(ib,iparm)*1.987D-3) ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0 #elif defined(FUNCT) ft(6)=quot #else ft(6)=1.0d0 #endif ! quotl=1.0d0 ! kfacl=1.0d0 ! do l=1,5 ! quotl1=quotl ! quotl=quotl*quot ! kfacl=kfacl*kfac ! fT(l)=kfacl/(kfacl-1.0d0+quotl) ! enddo ! else if (rescale_mode.eq.2) then ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) !#if defined(FUNCTH) ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3) ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0 !#elif defined(FUNCT) ! ft(6)=quot !#else ! ft(6)=1.0d0 !#endif ! quotl=1.0d0 ! do l=1,5 ! quotl=quotl*quot ! fT(l)=1.12692801104297249644d0/ & ! dlog(dexp(quotl)+dexp(-quotl)) ! enddo ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft ! else if (rescale_mode.eq.0) then ! do l=1,5 ! fT(l)=0.0d0 ! enddo ! else ! write (iout,*) & ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode ! call flush(iout) ! return 1 ! endif ! el end old rescale weihgts call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3)) #ifdef MPI do i=1,scount(me1) #else do i=1,ntot(islice) #endif evdw=enetb(1,i,iparm) ! evdw_t=enetb(21,i,iparm) evdw_t=enetb(20,i,iparm) #ifdef SCP14 ! evdw2_14=enetb(17,i,iparm) evdw2_14=enetb(18,i,iparm) evdw2=enetb(2,i,iparm)+evdw2_14 #else evdw2=enetb(2,i,iparm) evdw2_14=0.0d0 #endif #ifdef SPLITELE ees=enetb(3,i,iparm) evdw1=enetb(16,i,iparm) #else ees=enetb(3,i,iparm) evdw1=0.0d0 #endif ecorr=enetb(4,i,iparm) ecorr5=enetb(5,i,iparm) ecorr6=enetb(6,i,iparm) eel_loc=enetb(7,i,iparm) eello_turn3=enetb(8,i,iparm) eello_turn4=enetb(9,i,iparm) eello_turn6=enetb(10,i,iparm) ebe=enetb(11,i,iparm) escloc=enetb(12,i,iparm) etors=enetb(13,i,iparm) etors_d=enetb(14,i,iparm) ehpb=enetb(15,i,iparm) estr=enetb(17,i,iparm) ! estr=enetb(18,i,iparm) ! esccor=enetb(19,i,iparm) esccor=enetb(21,i,iparm) ! edihcnstr=enetb(20,i,iparm) edihcnstr=enetb(19,i,iparm) ecationcation=enetb(42,i,iparm) ecation_prot=enetb(41,i,iparm) !#ifdef SPLITELE ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & ! +wvdwpp*evdw1 & ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc & ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & ! +ft(2)*wturn3*eello_turn3 & ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & ! +wbond*estr !#else ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & ! +ft(1)*welec*(ees+evdw1) & ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc & ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & ! +ft(2)*wturn3*eello_turn3 & ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & ! +wbond*estr !#endif #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4 & +wturn3*eello_turn3 & +wturn6*eello_turn6+wel_loc*eel_loc & +edihcnstr+wtor_d*etors_d+wsccor*esccor & +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation #else etot=wsc*evdw+wscp*evdw2 & +welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4 & +wturn3*eello_turn3 & +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr & +wtor_d*etors_d+wsccor*esccor & +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation #endif #ifdef MPI Fdimless(i)= & beta_h(ib,iparm)*etot-entfac(i) potE(i,iparm)=etot #ifdef DEBUG write (iout,*) i,indstart(me)+i-1,ib,& 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),& -entfac(i),Fdimless(i) #endif #else Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i) potE(i,iparm)=etot #endif enddo ! i #ifdef MPI do i=1,scount(me1) Fdimless_(i)=Fdimless(i) enddo call MPI_Gatherv(Fdimless_(1),scount(me),& MPI_REAL,Fdimless(1),& scount_(0),idispl(0),MPI_REAL,Master,& WHAM_COMM, IERROR) #ifdef DEBUG call MPI_Gatherv(potE(1,iparm),scount_(me),& MPI_DOUBLE_PRECISION,potE(1,iparm),& scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,& WHAM_COMM, IERROR) call MPI_Gatherv(entfac(1),scount(me),& MPI_DOUBLE_PRECISION,entfac(1),& scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,& WHAM_COMM, IERROR) #endif if (me.eq.Master) then #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ntot(islice) write (iout,*) i,fdimless(i) enddo #endif #endif do i=1,ntot(islice) iperm(i)=i enddo call mysort1(ntot(islice),Fdimless,iperm) #ifdef DEBUG write (iout,*) "The FDIMLESS array after sorting" do i=1,ntot(islice) write (iout,*) i,iperm(i),fdimless(i) enddo #endif qfree=0.0d0 do i=1,ntot(islice) qfree=qfree+exp(-fdimless(i)+fdimless(1)) enddo ! write (iout,*) "qfree",qfree nlist=1 sumprob=0.0 do i=1,min0(ntot(islice),ensembles) sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree #ifdef DEBUG write (iout,*) i,ib,beta_h(ib,iparm),& 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),& potE(iperm(i),iparm),& -entfac(iperm(i)),fdimless(i),sumprob #endif if (sumprob.gt.0.99d0) goto 122 nlist=nlist+1 enddo 122 continue #ifdef MPI endif call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,& IERROR) call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,& IERROR) do i=1,nlist ii=iperm(i) iproc=0 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc)) iproc=iproc+1 enddo if (iproc.ge.nprocs) then write (iout,*) "Fatal error: processor out of range",iproc call flush(iout) if (bxfile) then close (ientout) else close (ientout,status="delete") endif return 1 endif ik=ii-indstart(iproc)+1 if (iproc.ne.Master) then if (me.eq.iproc) then #ifdef DEBUG write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,& " energy",potE(ik,iparm) #endif call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,& Master,i,WHAM_COMM,IERROR) else if (me.eq.Master) then call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,& WHAM_COMM,STATUS,IERROR) endif else if (me.eq.Master) then enepot(i)=potE(ik,iparm) endif enddo #else do i=1,nlist enepot(i)=potE(iperm(i),iparm) enddo #endif #ifdef MPI if (me.eq.Master) then #endif write(licz3,'(bz,i3.3)') iparm write(licz2,'(bz,i2.2)') islice if (temper.lt.100.0d0) then write(ctemper,'(f3.0)') temper else if (temper.lt.1000.0) then write (ctemper,'(f4.0)') temper else write (ctemper,'(f5.0)') temper endif if (nparmset.eq.1) then if (separate_parset) then write(licz4,'(bz,i3.3)') myparm pdbname=prefix(:ilen(prefix))//"_par"//licz4 else pdbname=prefix(:ilen(prefix)) endif else pdbname=prefix(:ilen(prefix))//"_parm_"//licz3 endif if (nslice.eq.1) then pdbname=pdbname(:ilen(pdbname))//"_T_"// & ctemper(:ilen(ctemper))//"pdb" else pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// & ctemper(:ilen(ctemper))//"pdb" endif open(ipdb,file=pdbname) do i=1,nlist read (ientout,rec=iperm(i)) & ((csingle(l,k),l=1,3),k=1,nres),& ((csingle(l,k+nres),l=1,3),k=nnt,nct),& nss,(ihpb(k),jhpb(k),k=1,nss),& eini,efree,rmsdev,iscor do j=1,2*nres do k=1,3 c(k,j)=csingle(k,j) enddo enddo call returnbox do j=1,2*nres do k=1,3 csingle(k,j)=c(k,j) enddo enddo eini=fdimless(i) call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev) enddo #ifdef MPI endif #endif enddo ! ib enddo ! iparm if (bxfile) then close(ientout) else close(ientout,status="delete") endif do i=0,nprocs scount(i)=scount_(i) enddo return end subroutine make_ensembles !-------------------------------------------------------------------------- subroutine mysort1(n, x, ipermut) ! implicit none integer :: i,j,imax,ipm,n real(kind=4) :: x(n) integer :: ipermut(n) real(kind=4) :: xtemp do i=1,n xtemp=x(i) imax=i do j=i+1,n if (x(j).lt.xtemp) then imax=j xtemp=x(j) endif enddo x(imax)=x(i) x(i)=xtemp ipm=ipermut(imax) ipermut(imax)=ipermut(i) ipermut(i)=ipm enddo return end subroutine mysort1 !-------------------------------------------------------------------------- subroutine alloc_enecalc_arrays(iparm) use control_data use geometry_data, only:maxlob integer :: iparm !--------------------------- ! COMMON.ENERGIES form wham_data ! common /energies/ allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm) allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc) allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc) allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm) ! ! allocate ENECALC arrays !--------------------------- ! COMMON.ALLPARM ! common /allparm/ allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm) allocate(vbldsc0_all(maxbondterm,ntyp,iparm),& aksc_all(maxbondterm,ntyp,iparm),& abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm) allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm) allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),& bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm) allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm) allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm) allocate(theta0_all(-ntyp:ntyp,iparm),& sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm) allocate(aa0thet_all(-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm)) !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm) allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm)) !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm) allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm)) allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm)) allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm)) allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm)) !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1, ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm) allocate(ffthet_all1(ndouble,ndouble,ntheterm3,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,iparm)) allocate(ggthet_all1(ndouble,ndouble,ntheterm3,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,iparm)) allocate(ffthet_all2(ndouble,ndouble,ntheterm3,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,iparm)) allocate(ggthet_all2(ndouble,ndouble,ntheterm3,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,iparm)) !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm) allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm) allocate(bsc_all(maxlob,ntyp,iparm)) !(maxlob,ntyp,max_parm) allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm) allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm) allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm) allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm)) allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm)) allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm) allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,& -ntortyp:ntortyp,2,iparm)) allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,& -ntortyp:ntortyp,2,iparm)) !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,& -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,& -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(b1_all(2,-nloctyp:nloctyp,iparm)) allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm) allocate(cc_all(2,2,-nloctyp:nloctyp,iparm)) allocate(dd_all(2,2,-nloctyp:nloctyp,iparm)) allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm) allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm)) allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm) allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm) allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),& ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm) allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm) allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),& augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),& sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),& chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm) allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm) allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),& akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),& v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm) allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm)) allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm) allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm) allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm)) !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,& -ntortyp:ntortyp,2,iparm)) allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,& -ntortyp:ntortyp,2,iparm)) !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm) allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm) allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm) allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),& ntheterm2_all(iparm),ntheterm3_all(nParmSet),& nsingle_all(iparm),& ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm) allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm) ! end subroutine alloc_enecalc_arrays !-------------------------------------------------------------------------- !-------------------------------------------------------------------------- end module ene_calc