subroutine w2x(nvarr,x,*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CONTROL" include "COMMON.INTERACT" include "COMMON.TORSION" include "COMMON.LOCAL" 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/ integer ind_shield /25/ double precision x(max_paropt),xchuj,epsi integer nvarr,i,ii,j,k,l,ll,iprot,iblock integer iti,itj c write (iout,*) "W2X: SPLIT_FOURIERTOR",SPLIT_FOURIERTOR, c & " tor_mode",tor_mode ii=0 do i=1,n_ene if (imask(i).eq.1 .and. i.ne.ind_shield) then ii=ii+1 x(ii)=ww(i) x_low(ii)=ww_low(i) x_up(ii)=ww_up(i) endif enddo #ifdef NEWCORR c 2/10/208 AL: insert the base torsional PMFs if (mod_fourier(nloctyp).gt.0) then if (torsion_pmf) then do iti=0,2 do itj=-2,2 ii=ii+1 x(ii)=e0(iti,itj) x_low(ii)=e0_low(iti,itj) x_up(ii)=e0_up(iti,itj) mask_e0(iti,itj)=ii enddo enddo else do iti=0,2 do itj=-2,2 mask_e0(iti,itj)=0 enddo enddo endif if (torsion_pmf.and.eello_pmf) then ii=ii+1 x(ii)=wello_PMF x_low(ii)=wello_PMF_low x_up(ii)=wello_PMF_up mask_ello_PMF=ii else mask_ello_PMF=0 endif if (turn_pmf.and.(torsion_pmf.or.eello_pmf)) then ii=ii+1 x(ii)=wturn_PMF x_low(ii)=wturn_PMF_low x_up(ii)=wturn_PMF_up mask_turn_PMF=ii else mask_turn_PMF=0 endif endif #endif n_weight_PMF=ii 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 c write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc, c & " npair_sc",npair_sc c 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 ntor_var=0 nbacktor_var=0 if (mod_tor) then c 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) c write (iout,'(2a4,2i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,3i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,3i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,3i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,2i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,2i4)') restyp(itor2typ(j)), c & 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) c write (iout,'(2a4,2i4)') restyp(itor2typ(j)), c & restyp(itor2typ(i)),k,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 c write (iout,*) "i",i," j",j," weitor",weitor(0,j,i,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) endif #else write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR" stop #endif endif ! tor_mode enddo enddo endif nbacktor_var=ntor_var c write (iout,*) "nbacktor_var",nbacktor_var if (mod_sccor) then c 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) c 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) c 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) c 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) c write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii enddo endif enddo enddo enddo endif c 7/8/17 AL Optimization of the bending potentials nang_var=0 if (mod_ang) then do i=0,nthetyp-1 if (mask_ang(i).eq.0) cycle do j=1,nbend_kcc_TB(i) ii=ii+1 nang_var=nang_var+1 x(ii)=v1bend_chyb(j,i) x_low(ii)=v1bend_low(j,i) x_up(ii)=v1bend_up(j,i) enddo enddo endif c write (iout,*) "nang_var",nang_var c Optimized cumulant coefficients 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 mask_bnew1(k,j,i)=ii if (mask_bnewtor1(k,j,i).eq.0) mask_bnewtor1(k,j,i)=-ii 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 mask_bnew2(k,j,i)=ii 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 mask_ccnew(k,j,i)=ii 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 mask_ddnew(k,j,i)=ii 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,3 if (mask_e0new(j,i).gt.0) then ii=ii+1 mask_e0new(j,i)=ii 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 mask_eenew(l,k,j,i)=ii 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 IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN do i=0,nloctyp-1 c write (iout,*) "type",i," mod_fouriertor",mod_fouriertor(i) if (mod_fouriertor(i).gt.0) then do j=1,2 do k=1,3 c write (iout,*) "mask_bnew2tor",k,j,mask_bnew1tor(k,j,i) if (mask_bnew1tor(k,j,i).gt.0) then ii=ii+1 mask_bnew1tor(k,j,i)=ii x(ii)=bnew1tor(k,j,i) x_low(ii)=bnew1tor_low(k,j,i) x_up(ii)=bnew1tor_up(k,j,i) endif enddo enddo do j=1,2 do k=1,3 c write (iout,*) "mask_bnew2tor",j,k,mask_bnew2tor(k,j,i) if (mask_bnew2tor(k,j,i).gt.0) then ii=ii+1 mask_bnew2tor(k,j,i)=ii x(ii)=bnew2tor(k,j,i) x_low(ii)=bnew2tor_low(k,j,i) x_up(ii)=bnew2tor_up(k,j,i) endif enddo enddo do j=1,2 do k=1,3 c write (iout,*) "mask_ccnewtor",k,j,mask_ccnewtor(k,j,i) if (mask_ccnewtor(k,j,i).gt.0) then ii=ii+1 mask_ccnewtor(k,j,i)=ii x(ii)=ccnewtor(k,j,i) x_low(ii)=ccnewtor_low(k,j,i) x_up(ii)=ccnewtor_up(k,j,i) endif enddo enddo do j=1,2 do k=1,3 c write (iout,*) "mask_ddnewtor",k,j,mask_ddnewtor(k,j,i) if (mask_ddnewtor(k,j,i).gt.0) then ii=ii+1 mask_ddnewtor(k,j,i)=ii x(ii)=ddnewtor(k,j,i) x_low(ii)=ddnewtor_low(k,j,i) x_up(ii)=ddnewtor_up(k,j,i) endif enddo enddo c write (iout,*) "type",i do j=1,3 c write (iout,*) "mask_e0newtor",j,mask_e0newtor(j,i) if (mask_e0newtor(j,i).gt.0) then ii=ii+1 mask_e0newtor(j,i)=ii x(ii)=e0newtor(j,i) x_low(ii)=e0newtor_low(j,i) x_up(ii)=e0newtor_up(j,i) c write (iout,*) j," e0newtorlow",e0newtor_low(j,i), c & " e0newtorup",e0newtor_up(j,i) endif enddo do j=1,2 do k=1,2 do l=1,2 c write (iout,*)"mask_eenewtor",l,k,j,mask_eenewtor(l,k,j,i) if (mask_eenewtor(l,k,j,i).gt.0) then ii=ii+1 mask_eenewtor(l,k,j,i)=ii x(ii)=eenewtor(l,k,j,i) x_low(ii)=eenewtor_low(l,k,j,i) x_up(ii)=eenewtor_up(l,k,j,i) c write (iout,*) j,k,l, c & " eenewtorlow",eenewtor_low(l,k,j,i), c & " eenewtorup",eenewtor_up(l,k,j,i) endif enddo enddo enddo endif enddo ELSE do i=0,nloctyp-1 do j=1,2 do k=1,3 mask_bnewtor1(k,j,i)=-mask_bnew1(k,j,i) enddo enddo do j=1,2 do k=1,3 mask_bnewtor2(k,j,i)=-mask_bnew2(k,j,i) enddo enddo do j=1,2 do k=1,3 mask_ccnewtor(k,j,i)=-mask_ccnew(k,j,i) enddo enddo do j=1,2 do k=1,3 mask_ddnewtor(k,j,i)=-mask_ddnew(k,j,i) enddo enddo do j=1,3 mask_e0newtor(j,i)=-mask_e0new(j,i) enddo do j=1,2 do k=1,2 do l=1,2 mask_eenewtor(l,k,j,i)=-mask_eenew(l,k,j,i) enddo enddo enddo enddo ENDIF #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 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 if (mask_sigma(5,i).gt.0) then ii=ii+1 x(ii)=rr0(i) x_low(ii)=sigma_low(5,i) x_up(ii)=sigma_up(5,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 Weight of the shielding term if (imask(ind_shield).eq.1) then ii=ii+1 x(ii)=ww(ind_shield) x_low(ii)=ww_low(ind_shield) x_up(ii)=ww_up(ind_shield) 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,3a15)') " x"," x_low"," x_up" do i=1,nvarr write (iout,'(i3,3e15.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.LOCAL" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.CLASSES" include "COMMON.WEIGHTDER" include "COMMON.SCCOR" include "COMMON.SHIELD" 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/ integer ind_shield /25/ lprint=.false. if (lprint) write (iout,*) "x2w: nvar",nvarr if (lprint) write (iout,*) "x",(x(i),i=1,nvarr) if (lprint) write(iout,*)"X2W: SPLIT_FOURIERTOR",SPLIT_FOURIERTOR, & " tor_mode",tor_mode ii=0 do i=1,n_ene if (imask(i).eq.1 .and. i.ne.ind_shield) then ii=ii+1 ww(i)=x(ii) c write (iout,*) i, wname(i)," ",ww(i) endif enddo #ifdef NEWCORR c 2/10/208 AL: insert the base torsional PMFs if (mod_fourier(nloctyp).gt.0) then do iti=0,2 do itj=-2,2 if (mask_e0(iti,itj).gt.0) then ii=ii+1 e0(iti,itj)=x(ii) endif enddo enddo if (mask_ello_PMF.gt.0) then ii=ii+1 wello_PMF=x(ii) endif if (mask_turn_PMF.gt.0) then ii=ii+1 wturn_PMF=x(ii) endif endif #endif 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 #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 endif #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 7/8/17 AL Optimization of the bending potentials if (mod_ang) then do i=0,nthetyp-1 if (mask_ang(i).eq.0) cycle do j=1,nbend_kcc_TB(i) ii=ii+1 v1bend_chyb(j,i)=x(ii) v1bend_chyb(j,-i)=x(ii) enddo 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,3 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 IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN do i=0,nloctyp-1 if (mod_fouriertor(i).gt.0) then do j=1,2 do k=1,3 if (mask_bnew1tor(k,j,i).gt.0) then ii=ii+1 bnew1tor(k,j,i)=x(ii) endif enddo enddo do j=1,2 do k=1,3 if (mask_bnew2tor(k,j,i).gt.0) then ii=ii+1 bnew2tor(k,j,i)=x(ii) endif enddo enddo do j=1,2 do k=1,3 if (mask_ccnewtor(k,j,i).gt.0) then ii=ii+1 ccnewtor(k,j,i)=x(ii) endif enddo enddo do j=1,2 do k=1,3 if (mask_ddnewtor(k,j,i).gt.0) then ii=ii+1 ddnewtor(k,j,i)=x(ii) endif enddo enddo do j=1,3 if (mask_e0newtor(j,i).gt.0) then ii=ii+1 e0newtor(j,i)=x(ii) endif enddo do j=1,2 do k=1,2 do l=1,2 if (mask_eenewtor(l,k,j,i).gt.0) then ii=ii+1 eenewtor(l,k,j,i)=x(ii) endif enddo enddo enddo endif enddo do i=1,nloctyp-1 do j=1,3 bnew1tor(j,1,-i)= bnew1tor(j,1,i) bnew1tor(j,2,-i)=-bnew1tor(j,2,i) bnew2tor(j,1,-i)= bnew2tor(j,1,i) bnew2tor(j,2,-i)=-bnew2tor(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 ccnewtor(j,1,-i)=ccnewtor(j,1,i) ccnewtor(j,2,-i)=-ccnewtor(j,2,i) ddnewtor(j,1,-i)=ddnewtor(j,1,i) ddnewtor(j,2,-i)=-ddnewtor(j,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)') & "Coefficients of the single-body torsional 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,(bnew1tor(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,(bnew2tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of C" 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 D" 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 E" 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 call flush(iout) endif ELSE do i=-nloctyp+1,nloctyp-1 do k=1,3 do j=1,2 bnew1tor(k,j,i)=bnew1(k,j,i) enddo enddo do k=1,3 do j=1,2 bnew2tor(k,j,i)=bnew2(k,j,i) enddo enddo do k=1,3 ccnewtor(k,1,i)=ccnew(k,1,i) ccnewtor(k,2,i)=ccnew(k,2,i) ddnewtor(k,1,i)=ddnew(k,1,i) ddnewtor(k,2,i)=ddnew(k,2,i) enddo do j=1,3 e0newtor(j,i)=e0new(j,i) enddo do j=1,2 do k=1,2 do l=1,2 eenewtor(l,k,j,i)=eenew(l,k,j,i) enddo enddo enddo enddo ENDIF ! SPLIT_FOURIERTOR if (tor_mode.eq.2) then c Recalculate new torsional potential coefficients do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 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) c v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) c & -ccnewtor(k,2,i)*ddnewtor(l,2,j)) c v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) c & +ccnewtor(k,1,i)*ddnewtor(l,2,j)) 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)) enddo enddo enddo enddo 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 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 if (mask_sigma(5,i).gt.0) then ii=ii+1 rr0(i)=x(ii) c write(iout,*) "alp",i,ii,x(ii) endif enddo endif c write (iout,*)"mod_side",mod_side," mod_side_ohter",mod_side_other 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) do i=1,ntyp chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0) enddo 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) then write (iout,'(/a)') 'One-body parameters:' if (ipot.eq.1) then write (iout,'(a,4x,a)') 'residue','sigma' write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp) else if (ipot.eq.2) then write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 ' write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i), & i=1,ntyp) else if (ipot.eq.3 .or. ipot.eq.4) then write (iout,'(a,4x,5a)') 'residue',' sigma ','s||/s_|_^2', & ' chip0 ',' chip ',' alph ' write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i), & chip0(i),chip(i),alp(i),i=1,ntyp) else if (ipot.eq.5) then write (iout,'(a,4x,6a)') 'residue',' sigma ',' r0 ', & 's||/s_|_^2',' chip0 ',' chip ',' alph ' write (iout,'(a3,6x,6f10.5)') (restyp(i),sigma0(i),rr0(i), & sigii(i),chip0(i),chip(i),alp(i),i=1,ntyp) endif write (iout,'(/a/10x,7a/72(1h-))') & 'Working parameters of the SC interactions:', & ' a ',' b ',' augm ',' sigma ',' r0 ', & ' chi1 ',' chi2 ' endif 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 Weight of the shielding term if (imask(ind_shield).eq.1) then ii=ii+1 ww(ind_shield)=x(ii) endif wshield=ww(ind_shield) 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