subroutine w2x(nvarr,x,*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CONTROL" include "COMMON.INTERACT" include "COMMON.TORSION" include "COMMON.WEIGHTS" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.NAMES" include "COMMON.VAR" include "COMMON.VMCPAR" include "COMMON.CLASSES" include "COMMON.WEIGHTDER" include "COMMON.FFIELD" include "COMMON.SCCOR" integer itor2typ(-2:2) /-20,-20,10,9,20/ double precision x(max_paropt),xchuj,epsi integer nvarr,i,ii,j,k,l,ll,iprot,iblock integer iti,itj ii=0 do i=1,n_ene if (imask(i).eq.1) then ii=ii+1 x(ii)=ww(i) x_low(ii)=ww_low(i) x_up(ii)=ww_up(i) endif enddo c c SC parameters c Change 9/21/07 AL Variables are epsilons and not the increments on epsilons c to enable control on their deviations from original values. c Single-body epsilons are now defined as epsilons of pairs of same type side c chains. c write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc, & " npair_sc",npair_sc call flush(iout) if (mod_side) then do i=1,ntyp do j=1,ntyp eps_orig(i,j)=eps(i,j) enddo enddo do i=1,nsingle_sc iti=ityp_ssc(i) ii=ii+1 epsi=eps(iti,iti) x(ii)=epsi c write (iout,*) "ii",ii," x",x(ii) x_low(ii)=epss_low(i) x_up(i)=epss_up(i) c if (epss_up(i).eq.0.0d0) then c x_up(ii)=dabs(epsi)*epss_low(i) c x_low(ii)=epsi-x_up(ii) c x_up(ii)=epsi+x_up(ii) c else c x_low(ii)=epss_low(i) c x_up(ii)=epss_up(i) c endif enddo do i=1,npair_sc ii=ii+1 iti=ityp_psc(1,i) itj=ityp_psc(2,i) if (iti.eq.itj) then do j=1,nsingle_sc if (ityp_ssc(j).eq.iti) then write (iout,*) "Error - pair ",restyp(iti),"-", & restyp(iti), & " specified in variable epsilon list but type defined", & " in single-type list." return1 endif enddo endif epsi=eps(iti,itj) x(ii)=epsi x_low(ii)=epsp_low(i) x_up(ii)=epsp_up(i) c if (epsp_up(i).eq.0.0d0) then c x_up(ii)=dabs(epsi)*epsp_low(i) c x_low(ii)=epsi-x_up(ii) c x_up(ii)=epsi+x_up(ii) c else c x_low(ii)=epsp_low(i) c x_up(ii)=epsp_up(i) c endif c write (iout,*) "ii",ii," x",x(ii) enddo endif if (mod_fourier(nloctyp).gt.0) then #ifdef NEWCORR do i=0,nloctyp-1 if (mod_fourier(i).gt.0) then do j=1,2 do k=1,3 if (mask_bnew1(k,j,i).gt.0) then ii=ii+1 x(ii)=bnew1(k,j,i) x_low(ii)=bnew1_low(k,j,i) x_up(ii)=bnew1_up(k,j,i) endif enddo enddo do j=1,2 do k=1,3 if (mask_bnew2(k,j,i).gt.0) then ii=ii+1 x(ii)=bnew2(k,j,i) x_low(ii)=bnew2_low(k,j,i) x_up(ii)=bnew2_up(k,j,i) endif enddo enddo do j=1,2 do k=1,3 if (mask_ccnew(k,j,i).gt.0) then ii=ii+1 x(ii)=ccnew(k,j,i) x_low(ii)=ccnew_low(k,j,i) x_up(ii)=ccnew_up(k,j,i) endif enddo enddo do j=1,2 do k=1,3 if (mask_ddnew(k,j,i).gt.0) then ii=ii+1 x(ii)=ddnew(k,j,i) x_low(ii)=ddnew_low(k,j,i) x_up(ii)=ddnew_up(k,j,i) endif enddo enddo do j=1,2 if (mask_e0new(j,i).gt.0) then ii=ii+1 x(ii)=e0new(j,i) x_low(ii)=e0new_low(j,i) x_up(ii)=e0new_up(j,i) endif enddo do j=1,2 do k=1,2 do l=1,2 if (mask_eenew(l,k,j,i).gt.0) then ii=ii+1 x(ii)=eenew(l,k,j,i) x_low(ii)=eenew_low(l,k,j,i) x_up(ii)=eenew_up(l,k,j,i) endif enddo enddo enddo endif enddo #else do i=1,nloctyp if (mod_fourier(i).gt.0) then do j=2,13 if (mask_fourier(j,i).gt.0) then ii=ii+1 x(ii)=b(j,i) x_low(ii)=b_low(j,i) x_up(ii)=b_up(j,i) endif enddo endif enddo #endif endif ntor_var=0 nbacktor_var=0 if (mod_tor) then write (iout,*) "Indices of optimizable torsionals" do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 c write (iout,*) "mask",i,j,mask_tor(j,i) if (tor_mode.eq.0) then do iblock=1,2 if (mask_tor(0,j,i,iblock).eq.1) then ii=ii+1 ntor_var=ntor_var+1 x(ii)=weitor(0,j,i,iblock) x_low(ii)=weitor_low(0,j,i,iblock) x_up(ii)=weitor_up(0,j,i,iblock) write (iout,'(2a4,2i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),iblock,ii else if (mask_tor(0,j,i,iblock).eq.2) then do k=1,nterm(j,i,iblock) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v1(k,j,i,iblock) x_low(ii)=weitor_low(0,j,i,iblock) x_up(ii)=weitor_up(0,j,i,iblock) write (iout,'(2a4,3i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),iblock,k,ii enddo else if (mask_tor(0,j,i,iblock).gt.2) then do k=1,nterm(j,i,iblock) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v1(k,j,i,iblock) x_low(ii)=weitor_low(0,j,i,iblock) x_up(ii)=weitor_up(0,j,i,iblock) write (iout,'(2a4,3i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),iblock,k,ii enddo do k=1,nterm(j,i,iblock) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v2(k,j,i,iblock) x_low(ii)=weitor_low(0,j,i,iblock) x_up(ii)=weitor_up(0,j,i,iblock) write (iout,'(2a4,3i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),ii,k,ii enddo endif enddo ! iblock else if (tor_mode.eq.1) then if (mask_tor(0,j,i,1).eq.1) then ii=ii+1 ntor_var=ntor_var+1 x(ii)=weitor(0,j,i,1) x_low(ii)=weitor_low(0,j,i,1) x_up(ii)=weitor_up(0,j,i,1) write (iout,'(2a4,i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),ii else if (mask_tor(0,j,i,1).eq.2) then do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v1_kcc(ll,l,k,j,i) x_low(ii)=weitor_low(0,j,i,1) x_up(ii)=weitor_up(0,j,i,1) write (iout,'(2a4,2i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),k,ii enddo ! ll enddo ! l enddo else if (mask_tor(0,j,i,1).gt.2) then do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v1_kcc(ll,l,k,j,i) x_low(ii)=weitor_low(0,j,i,1) x_up(ii)=weitor_up(0,j,i,1) write (iout,'(2a4,2i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),k,ii enddo ! ll enddo ! l enddo do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v2_kcc(ll,l,k,j,i) x_low(ii)=weitor_low(0,j,i,1) x_up(ii)=weitor_up(0,j,i,1) write (iout,'(2a4,2i4)') restyp(itor2typ(j)), & restyp(itor2typ(i)),k,ii enddo! ll enddo ! l enddo endif #ifndef NEWCORR else write (iout,*) "ERROR! TOR_MODE>1 ALLOWED ONLY WITH NEWCORR" #endif endif ! tor_mode enddo enddo endif nbacktor_var=ntor_var write (iout,*) "nbacktor_var",nbacktor_var if (mod_sccor) then write (iout,*) "Indices of optimizable SCCOR torsionals" do i=-nsccortyp,nsccortyp do j=-nsccortyp,nsccortyp do k=1,3 c write (iout,*) "mask",i,j,mask_tor(j,i) if (mask_tor(k,j,i,1).eq.1) then ii=ii+1 ntor_var=ntor_var+1 x(ii)=weitor(k,j,i,1) x_low(ii)=weitor_low(k,j,i,1) x_up(ii)=weitor_up(k,j,i,1) write (iout,'(i4,2a4,i4)') k,restyp(j),restyp(i),ii else if (mask_tor(k,j,i,1).eq.2) then do l=1,nterm_sccor(j,i) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v1sccor(l,k,j,i) x_low(ii)=weitor_low(k,j,i,1) x_up(ii)=weitor_up(k,j,i,1) write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii enddo else if (mask_tor(k,j,i,1).gt.2) then do l=1,nterm_sccor(j,i) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v1sccor(l,k,j,i) x_low(ii)=weitor_low(k,j,i,1) x_up(ii)=weitor_up(k,j,i,1) write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii enddo do l=1,nterm_sccor(j,i) ii=ii+1 ntor_var=ntor_var+1 x(ii)=v2sccor(l,k,j,i) x_low(ii)=weitor_low(k,j,i,1) x_up(ii)=weitor_up(k,j,i,1) write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii enddo endif enddo enddo enddo endif c Added SC sigmas and anisotropies if (mod_side_other) then do i=1,ntyp if (mask_sigma(1,i).gt.0) then ii=ii+1 x(ii)=sigma0(i) x_low(ii)=sigma_low(1,i) x_up(ii)=sigma_up(1,i) c write (iout,*) "sig0",i,ii,x(ii), x_low(ii),x_up(ii) endif if (mask_sigma(2,i).gt.0) then ii=ii+1 x(ii)=sigii(i) x_low(ii)=sigma_low(2,i) x_up(ii)=sigma_up(2,i) c write (iout,*) "sigi",i,ii,x(ii), x_low(ii),x_up(ii) endif if (mask_sigma(3,i).gt.0) then ii=ii+1 x(ii)=chip0(i) x_low(ii)=sigma_low(3,i) x_up(ii)=sigma_up(3,i) c write (iout,*) "chi",i,ii,x(ii), x_low(ii),x_up(ii) endif if (mask_sigma(4,i).gt.0) then ii=ii+1 x(ii)=alp(i) x_low(ii)=sigma_low(4,i) x_up(ii)=sigma_up(4,i) c write (iout,*) "alp",i,ii,x(ii), x_low(ii),x_up(ii) endif enddo endif if (mod_elec) then do i=1,2 do j=1,i if (mask_elec(i,j,1).gt.0) then ii=ii+1 x(ii)=epp(i,j) x_low(ii)=epp_low(i,j) x_up(ii)=epp_up(i,j) endif if (mask_elec(i,j,2).gt.0) then ii=ii+1 x(ii)=rpp(i,j) x_low(ii)=rpp_low(i,j) x_up(ii)=rpp_up(i,j) endif if (mask_elec(i,j,3).gt.0) then ii=ii+1 x(ii)=elpp6(i,j) x_low(ii)=elpp6_low(i,j) x_up(ii)=elpp6_up(i,j) endif if (mask_elec(i,j,4).gt.0) then ii=ii+1 x(ii)=elpp3(i,j) x_low(ii)=elpp3_low(i,j) x_up(ii)=elpp3_up(i,j) endif enddo enddo endif C C Define the SC-p interaction constants C if (mod_scp) then if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+ & mask_scp(0,2,2) .gt. 0) then do j=1,2 if (mask_scp(0,j,1).gt.0) then ii=ii+1 x(ii)=eps_scp(1,j) x_low(ii)=epscp_low(0,j) x_up(ii)=epscp_up(0,j) endif if (mask_scp(0,j,2).gt.0) then ii=ii+1 x(ii)=rscp(1,j) x_low(ii)=rscp_low(0,j) x_up(ii)=rscp_up(0,j) endif enddo else do i=1,20 do j=1,2 if (mask_scp(i,j,1).gt.0) then ii=ii+1 x(ii)=eps_scp(i,j) x_low(ii)=epscp_low(i,j) x_up(ii)=epscp_up(i,j) endif if (mask_scp(i,j,2).gt.0) then ii=ii+1 x(ii)=rscp(i,j) x_low(ii)=rscp_low(i,j) x_up(ii)=rscp_up(i,j) endif enddo enddo endif endif c Compute boundaries, if defined relative to starting values do i=1,ii write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i) if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then xchuj=x_low(i) x_low(i)=x(i)-xchuj x_up(i)=x(i)+xchuj else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then xchuj=x_low(i) if (x(i).gt.0.0d0) then x_low(i)=x(i)*(1.0d0-xchuj) x_up(i)=x(i)*(1.0d0+xchuj) else x_low(i)=x(i)*(1.0d0+xchuj) x_up(i)=x(i)*(1.0d0-xchuj) endif endif enddo c Base of heat capacity do iprot=1,nprot if (nbeta(2,iprot).gt.0) then ii=ii+1 x(ii)=heatbase(iprot) x_low(ii)=-10.0d0 x_up(ii)=10.0d0 endif enddo nvarr=ii write (iout,*) "Number of optimizable parameters:",nvarr write (iout,*) write (iout,*) "Initial variables and boundaries:" write (iout,'(3x,3a10)') "x","x_low","x_up" do i=1,nvarr write (iout,'(i3,3f10.5)') i,x(i),x_low(i),x_up(i) enddo call flush(iout) return end c------------------------------------------------------------------------------ subroutine x2w(nvarr,x) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CONTROL" include "COMMON.NAMES" include "COMMON.WEIGHTS" include "COMMON.INTERACT" include "COMMON.ENERGIES" include "COMMON.FFIELD" include "COMMON.TORSION" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.CLASSES" include "COMMON.WEIGHTDER" include "COMMON.SCCOR" logical lprint integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock double precision rri,v0ij,si double precision x(nvarr) double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2, & ratsig1,ratsig2,rsum_max,dwa16,r_augm double precision eps_temp(ntyp,ntyp) double precision ftune_eps,ftune_epsprim logical in_sc_pair_list external in_sc_pair_list integer itor2typ(3) /10,9,20/ lprint=.false. if (lprint) write (iout,*) "x2w: nvar",nvarr if (lprint) write (iout,*) "x",(x(i),i=1,nvarr) ii=0 do i=1,n_ene if (imask(i).eq.1) then ii=ii+1 ww(i)=x(ii) endif enddo wsc=ww(1) wscp=ww(2) welec=ww(3) wcorr=ww(4) wcorr5=ww(5) wcorr6=ww(6) wel_loc=ww(7) wturn3=ww(8) wturn4=ww(9) wturn6=ww(10) wang=ww(11) wscloc=ww(12) wtor=ww(13) wtor_d=ww(14) wstrain=ww(15) wvdwpp=ww(16) wbond=ww(17) wsccor=ww(19) #ifdef SCP14 scal14=ww(18)/ww(2) #endif do i=1,n_ene weights(i)=ww(i) enddo if (lprint) then write (iout,*) "Weights:" do i=1,n_ene write (iout,'(a,f10.5)') wname(i),ww(i) enddo endif C SC parameters if (mod_side) then do i=1,nsingle_sc ii=ii+1 iti=ityp_ssc(i) eps(iti,iti)=x(ii) do j=1,iti-1 if (.not. in_sc_pair_list(iti,j)) then eps(iti,j)=eps_orig(iti,j) & +0.5d0*(x(ii)-eps_orig(iti,iti)) eps(j,iti)=eps(iti,j) endif enddo enddo do i=1,npair_sc ii=ii+1 iti=ityp_psc(1,i) itj=ityp_psc(2,i) eps(iti,itj)=x(ii) eps(itj,iti)=x(ii) enddo endif if (mod_fourier(nloctyp).gt.0) then #ifdef NEWCORR do i=0,nloctyp-1 if (mod_fourier(i).gt.0) then do j=1,2 do k=1,3 if (mask_bnew1(k,j,i).gt.0) then ii=ii+1 bnew1(k,j,i)=x(ii) endif enddo enddo do j=1,2 do k=1,3 if (mask_bnew2(k,j,i).gt.0) then ii=ii+1 bnew2(k,j,i)=x(ii) endif enddo enddo do j=1,2 do k=1,3 if (mask_ccnew(k,j,i).gt.0) then ii=ii+1 ccnew(k,j,i)=x(ii) endif enddo enddo do j=1,2 do k=1,3 if (mask_ddnew(k,j,i).gt.0) then ii=ii+1 ddnew(k,j,i)=x(ii) endif enddo enddo do j=1,2 if (mask_e0new(j,i).gt.0) then ii=ii+1 e0new(j,i)=x(ii) endif enddo do j=1,2 do k=1,2 do l=1,2 if (mask_eenew(l,k,j,i).gt.0) then ii=ii+1 eenew(l,k,j,i)=x(ii) endif enddo enddo enddo endif enddo do i=1,nloctyp-1 do j=1,3 bnew1(j,1,-i)= bnew1(j,1,i) bnew1(j,2,-i)=-bnew1(j,2,i) bnew2(j,1,-i)= bnew2(j,1,i) bnew2(j,2,-i)=-bnew2(j,2,i) enddo do j=1,3 c ccnew(j,1,i)=ccnew(j,1,i)/2 c ccnew(j,2,i)=ccnew(j,2,i)/2 c ddnew(j,1,i)=ddnew(j,1,i)/2 c ddnew(j,2,i)=ddnew(j,2,i)/2 ccnew(j,1,-i)=ccnew(j,1,i) ccnew(j,2,-i)=-ccnew(j,2,i) ddnew(j,1,-i)=ddnew(j,1,i) ddnew(j,2,-i)=-ddnew(j,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 call flush(iout) endif #else do i=0,nloctyp-1 if (mod_fourier(i).gt.0) then do j=2,13 if (mask_fourier(j,i).gt.0) then ii=ii+1 b(j,i)=x(ii) endif enddo if (lprint) then write (iout,*) "Fourier coefficients of cumulants" write (iout,*) 'Type',i write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=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,i) c B1(2,i) = b(5,i) c B1tilde(1,i) = b(3,i) c B1tilde(2,i) =-b(5,i) c B2(1,i) = b(2,i) c B2(2,i) = b(4,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) 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 Ctilde(1,1,i)=b(7,i) c Ctilde(1,2,i)=b(9,i) c Ctilde(2,1,i)=-b(9,i) c Ctilde(2,2,i)=b(7,i) DDold(1,1,i)= b(6,i) DDold(2,2,i)=-b(6,i) DDold(2,1,i)= b(8,i) DDold(1,2,i)= b(8,i) DDold(1,1,-i)= b(6,i) DDold(2,2,-i)=-b(6,i) DDold(2,1,-i)=-b(8,i) DDold(1,2,-i)=-b(8,i) c Dtilde(1,1,i)=b(6,i) c Dtilde(1,2,i)=b(8,i) c Dtilde(2,1,i)=-b(8,i) c Dtilde(2,2,i)=b(6,i) EEold(1,1,i)= b(10,i)+b(11,i) EEold(2,2,i)=-b(10,i)+b(11,i) EEold(2,1,i)= b(12,i)-b(13,i) EEold(1,2,i)= b(12,i)+b(13,i) EEold(1,1,i)= b(10,i)+b(11,i) EEold(2,2,i)=-b(10,i)+b(11,i) EEold(2,1,i)= b(12,i)-b(13,i) EEold(1,2,i)= b(12,i)+b(13,i) 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) endif 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 endif #ifdef NEWCORR if (mod_tor .or. mod_fourier(nloctyp).gt.0) then c write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier", c & mod_fourier(nloctyp) #else if (mod_tor) then #endif do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 if (tor_mode.eq.0) then do iblock=1,2 if (mask_tor(0,j,i,iblock).eq.1) then ii=ii+1 weitor(0,j,i,iblock)=x(ii) else if (mask_tor(0,j,i,iblock).eq.2) then do k=1,nterm(j,i,iblock) ii=ii+1 v1(k,j,i,iblock)=x(ii) enddo else if (mask_tor(0,j,i,iblock).gt.2) then do k=1,nterm(j,i,iblock) ii=ii+1 v1(k,j,i,iblock)=x(ii) enddo do k=1,nterm(j,i,iblock) ii=ii+1 v2(k,j,i,iblock)=x(ii) enddo endif v0ij=0.0d0 si=-1.0d0 do k=1,nterm(i,j,iblock) v1(k,-i,-j,iblock)=v1(k,i,j,iblock) v2(k,-i,-j,iblock)=-v2(k,i,j,iblock) v0ij=v0ij+si*v1(k,i,j,iblock) si=-si enddo do k=1,nlor(i,j,iblock) v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2) enddo v0(i,j,iblock)=v0ij v0(-i,-j,iblock)=v0ij enddo else if (tor_mode.eq.1) then if (mask_tor(0,j,i,1).eq.1) then ii=ii+1 weitor(0,j,i,1)=x(ii) else if (mask_tor(0,j,i,1).eq.2) then do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) ii=ii+1 v1_kcc(ll,l,k,j,i)=x(ii) enddo ! ll enddo ! l enddo else if (mask_tor(0,j,i,1).gt.2) then do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) ii=ii+1 v1_kcc(ll,l,k,j,i)=x(ii) enddo ! ll enddo ! l enddo do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) ii=ii+1 v2_kcc(ll,l,k,j,i)=x(ii) enddo! ll enddo ! l enddo endif else #ifdef NEWCORR c Compute torsional coefficients from local energy surface expansion coeffs if (mask_tor(0,j,i,1).eq.1) then ii=ii+1 weitor(0,j,i,1)=x(ii) c else if (mask_tor(0,j,i,1).eq.2) then else do k=1,nterm_kcc_Tb(j,i) do l=1,nterm_kcc_Tb(j,i) v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j) & +bnew1(k,2,i)*bnew2(l,2,j) v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j) & +bnew1(k,2,i)*bnew2(l,1,j) enddo enddo do k=1,nterm_kcc_Tb(j,i) do l=1,nterm_kcc_Tb(j,i) v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j) & -ccnew(k,2,i)*ddnew(l,2,j)) v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j) & +ccnew(k,1,i)*ddnew(l,2,j)) enddo enddo endif #else else write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR" stop #endif endif enddo enddo endif if (mod_sccor) then do i=-nsccortyp,nsccortyp do j=-nsccortyp,nsccortyp do k=1,3 if (mask_tor(k,j,i,1).eq.1) then ii=ii+1 weitor(k,j,i,1)=x(ii) else if (mask_tor(k,j,i,1).eq.2) then do l=1,nterm(j,i,1) ii=ii+1 v1sccor(l,k,j,i)=x(ii) enddo else if (mask_tor(k,j,i,1).gt.2) then do l=1,nterm_sccor(j,i) ii=ii+1 v1sccor(l,k,j,i)=x(ii) enddo do l=1,nterm_sccor(j,i) ii=ii+1 v2sccor(l,k,j,i)=x(ii) enddo endif enddo enddo enddo endif c SC sigmas and anisotropies as parameters if (mod_side_other) then c do i=1,ntyp c write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4) c enddo do i=1,ntyp if (mask_sigma(1,i).gt.0) then ii=ii+1 sigma0(i)=x(ii) c write(iout,*) "sig0",i,ii,x(ii) endif if (mask_sigma(2,i).gt.0) then ii=ii+1 sigii(i)=x(ii) c write(iout,*) "sigi",i,ii,x(ii) endif if (mask_sigma(3,i).gt.0) then ii=ii+1 chip0(i)=x(ii) c write(iout,*) "chip",i,ii,x(ii) endif if (mask_sigma(4,i).gt.0) then ii=ii+1 alp(i)=x(ii) c write(iout,*) "alp",i,ii,x(ii) endif enddo endif if (mod_side.or.mod_side_other) then if (lprint) write(iout,'(/3a,2i3)') 'Potential is ', & potname(ipot),', exponents are ',expon,2*expon C----------------------------------------------------------------------- C Calculate the "working" parameters of SC interactions. dwa16=2.0d0**(1.0d0/6.0d0) if (ipot.eq.4) then do i=1,ntyp chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0) enddo endif do i=1,ntyp do j=i,ntyp sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2) sigma(j,i)=sigma(i,j) rs0(i,j)=dwa16*sigma(i,j) rs0(j,i)=rs0(i,j) enddo enddo if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) if (lprint) write (iout,'(/a/10x,7a/72(1h-))') & 'Working parameters of the SC interactions:', & ' a ',' b ',' augm ',' sigma ',' r0 ', & ' chi1 ',' chi2 ' do i=1,ntyp do j=i,ntyp epsij=eps(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) else rrij=rr0(i)+rr0(j) endif r0(i,j)=rrij r0(j,i)=rrij rrij=rrij**expon epsij=eps(i,j) c sigeps=dsign(1.0D0,epsij) aa(i,j)=ftune_eps(epsij)*rrij*rrij bb(i,j)=-epsij*rrij aa(j,i)=aa(i,j) bb(j,i)=bb(i,j) if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 sigii1=sigii(i) sigii2=sigii(j) ratsig1=sigt2sq/sigt1sq ratsig2=1.0D0/ratsig1 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1) if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2) rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq) else rsum_max=sigma(i,j) endif sigmaii(i,j)=rsum_max sigmaii(j,i)=rsum_max if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) & then r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij augm(i,j)=epsij*r_augm**(2*expon) c augm(i,j)=0.5D0**(2*expon)*aa(i,j) augm(j,i)=augm(i,j) else augm(i,j)=0.0D0 augm(j,i)=0.0D0 endif if (lprint) then write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j), & sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo enddo endif if (mod_elec) then do i=1,2 do j=1,i if (mask_elec(i,j,1).gt.0) then ii=ii+1 epp(i,j)=x(ii) if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0 epp(j,i)=x(ii) endif if (mask_elec(i,j,2).gt.0) then ii=ii+1 rpp(i,j)=dabs(x(ii)) rpp(j,i)=dabs(x(ii)) endif if (mask_elec(i,j,3).gt.0) then ii=ii+1 elpp6(i,j)=x(ii) if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0 elpp6(j,i)=elpp6(i,j) endif if (mask_elec(i,j,4).gt.0) then ii=ii+1 elpp3(i,j)=x(ii) elpp3(j,i)=x(ii) endif rri=rpp(i,j)**6 app (i,j)=epp(i,j)*rri*rri app (j,i)=epp(j,i)*rri*rri bpp (i,j)=-2.0D0*epp(i,j)*rri bpp (j,i)=-2.0D0*epp(i,j)*rri ael6(i,j)=elpp6(i,j)*4.2D0**6 ael6(j,i)=elpp6(i,j)*4.2D0**6 ael3(i,j)=elpp3(i,j)*4.2D0**3 ael3(j,i)=elpp3(i,j)*4.2D0**3 enddo enddo if (lprint) then write (iout,'(/a)') 'Electrostatic interaction constants:' write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') & 'IT','JT','APP','BPP','AEL6','AEL3' do i=1,2 do j=1,2 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j), & ael6(i,j),ael3(i,j) enddo enddo endif endif C C Define the SC-p interaction constants C if (mod_scp) then if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+ & mask_scp(0,2,2) .gt. 0) then do j=1,2 if (mask_scp(0,j,1).gt.0) then ii=ii+1 do i=1,20 eps_scp(i,j)=x(ii) enddo endif if (mask_scp(0,j,2).gt.0) then ii=ii+1 do i=1,20 rscp(i,j)=x(ii) enddo endif enddo else do i=1,20 do j=1,2 if (mask_scp(i,j,1).gt.0) then ii=ii+1 eps_scp(i,j)=x(ii) endif if (mask_scp(i,j,2).gt.0) then ii=ii+1 rscp(i,j)=x(ii) endif enddo enddo endif do i=1,20 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6 enddo if (lprint) then write (iout,*) "Parameters of SC-p interactions:" do i=1,20 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1), & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2) enddo endif endif c write (iout,*) "x2w: nvar",nvarr c Base of heat capacity do iprot=1,nprot if (nbeta(2,iprot).gt.0) then ii=ii+1 heatbase(iprot)=x(ii) endif enddo if (ii.ne.nvarr) then write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)" write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii stop "CHUJ ABSOLUTNY!!!!" endif return end c--------------------------------------------------------------------------- double precision function ftune_eps(x) implicit none double precision x,y double precision delta /1.0d-5/ if (dabs(x).lt.delta) then y=x/delta ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8 else ftune_eps=dabs(x) endif return end c--------------------------------------------------------------------------- double precision function ftune_epsprim(x) implicit none double precision x,y double precision delta /1.0d-5/ if (dabs(x).lt.delta) then y=x/delta ftune_epsprim=1.5d0*y-0.5*y**3 else ftune_epsprim=dsign(1.0d0,x) endif return end c--------------------------------------------------------------------------- double precision function ftune_epsbis(x) implicit none double precision x,y double precision delta /1.0d-5/ if (dabs(x).lt.delta) then y=x/delta ftune_epsbis=(1.5d0-1.5d0*y**3)/delta else ftune_epsbis=0.0d0 endif return end c--------------------------------------------------------------------------- logical function in_sc_pair_list(iti,itj) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.WEIGHTS" integer i,iti,itj do i=1,npair_sc if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or. & ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then in_sc_pair_list=.true. return endif enddo in_sc_pair_list=.false. return end