subroutine minimize_init(n,iv,v,d) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" integer liv,lv parameter (liv=60,lv=(77+max_paropt*(max_paropt+17)/2)) include "COMMON.OPTIM" include "COMMON.MINPAR" include "COMMON.IOUNITS" include "COMMON.TIME1" integer n,iv(liv) double precision d(max_paropt),v(1:lv) integer i llocal=.false. call deflt(2,iv,liv,lv,v) * 12 means fresh start, dont call deflt iv(1)=12 * max num of fun calls iv(17)=maxfun * max num of iterations iv(18)=maxmin * controls output iv(19)=1 * selects output unit iv(21)=out_minim * 1 means to print out result iv(22)=print_fin * 1 means to print out summary stats iv(23)=print_stat * 1 means to print initial x and d iv(24)=print_ini * min val for v(radfac) default is 0.1 v(24)=0.1D0 * max val for v(radfac) default is 4.0 v(25)=2.0D0 c v(25)=4.0D0 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) * the sumsl default is 0.1 v(26)=0.1D0 * false conv if (act fnctn decrease) .lt. v(34) * the sumsl default is 100*machep v(34)=v(34)/100.0D0 * absolute convergence if (tolf.eq.0.0D0) tolf=1.0D-6 v(31)=tolf * relative convergence if (rtolf.eq.0.0D0) rtolf=1.0D-6 v(32)=rtolf * controls initial step size v(35)=1.0D-2 * large vals of d correspond to small components of step do i=1,n d(i)=1.0D0 enddo return end c------------------------------------------------------------------------------ subroutine minimize(n,x,f) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" integer liv,lv parameter (liv=60,lv=(77+max_paropt*(max_paropt+17)/2)) ********************************************************************* * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for * * the calling subprogram. * * when d(i)=1.0, then v(35) is the length of the initial step, * * calculated in the usual pythagorean way. * * absolute convergence occurs when the function is within v(31) of * * zero. unless you know the minimum value in advance, abs convg * * is probably not useful. * * relative convergence is when the model predicts that the function * * will decrease by less than v(32)*abs(fun). * ********************************************************************* #ifdef MPI include "mpif.h" integer IERROR include "COMMON.MPI" integer status(MPI_STATUS_SIZE) #endif include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.OPTIM" include "COMMON.CLASSES" include "COMMON.COMPAR" include "COMMON.MINPAR" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.TIME1" include "COMMON.ENERGIES" include "COMMON.PROTNAME" include "COMMON.THERMAL" integer nn common /cc/ nn integer ilen external ilen integer iv(liv) double precision x(max_paropt),d(max_paropt),v(1:lv),g(max_paropt) real p(max_paropt+1,max_paropt),y(max_paropt+1),pb(max_paropt), & pom(max_paropt),yb,y1,y2,pomi,ftol,temptr real funk external funk external targetfun,targetgrad,fdum integer uiparm,i,ii,j,k,nf,nfun,ngrad,iretcode,ianneal,maxanneal double precision urparm,penalty double precision f,f0,obf,obg(max_paropt),tcpu double precision Evarmax(maxprot),aux integer iscan,iprot,ibatch,niter integer n logical lprn lprn=.true. write (iout,*) "Minimizing with ",MINIMIZER IF (MINIMIZER.EQ."AMEBSA") THEN maxanneal=5 nn=n do i=1,n do j=1,n+1 p(j,i)=x(i) enddo enddo do j=1,n pom(j)=p(1,j) enddo y(1)=funk(pom) do i=2,n+1 do j=1,n pom(j)=p(i,j) enddo pom(i-1)=p(i,i-1)-0.2 pomi=pom(i-1) if (pom(i-1).lt.x_low(i-1)) pom(i-1)=x_low(i-1) y1=funk(pom) pom(i-1)=p(i,i-1)+0.2 if (pom(i-1).gt.x_up(i-1)) pom(i-1)=x_up(i-1) y2=funk(pom) if (y2.lt.y1) then p(i,i-1)=pom(i-1) y(i)=y2 else p(i,i-1)=pomi y(i)=y1 endif enddo yb=1.0e20 write (iout,*) "p and y" do i=1,n+1 write (iout,'(20e15.5)') (p(i,j),j=1,n),y(i) enddo ftol=rtolf do ianneal=1,maxanneal iretcode=maxmin temptr=5.0/ianneal write (iout,*) "ianneal",ianneal," temptr",temptr call amebsa(p,y,max_paropt+1,max_paropt,n,pb,yb,ftol,funk, & iretcode,temptr) write (iout,*) "Optimized P (PB)" do i=1,n write (iout,'(i5,f10.5)') i,pb(i) enddo write (iout,*) "Optimal function value",yb write(iout,*) "AMEBSA return code:",iretcode enddo do i=1,n x(i)=pb(i) enddo call targetfun(n,x,nf,f,uiparm,fdum) write (iout,*) write(iout,*) "Final function value:",f,f write(iout,*) "Final value of the object function:",obf write(iout,*) "Number of target function evaluations:",nfun write (iout,*) call x2w(n,x) do i=1, n xm(i)=x(i) enddo open (88,file=restartname,status="unknown") do i=1,n write (88,*) x(i) enddo close(88) ELSE IF (MINIMIZER.EQ."SUMSL") THEN nf=0 nfun=0 ngrad=0 niter=0 call minimize_init(n,iv,v,d) #ifdef CHECKGRAD write (iout,*) "Initial checkgrad:" call checkgrad(n,x) c write (iout,*) "Second checkgrad" c call checkgrad(n,x) #endif write (iout,*) "Calling SUMSL" iscan=0 call targetfun(n,x,nf,f,uiparm,fdum) write (iout,*) "Initial function value",f f0=f f=-1.0d20 do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 ) f0=f call scan(n,x,f,.false.) iscan=iscan+1 enddo PREVTIM=tcpu() f0=f f=-1.0d20 if (maxscan.gt.0) then iscan=0 do while ( iscan.lt.maxscan .and. (f0-f)/(f0+1.0d0).gt.0.2d0 ) 101 continue cutoffviol=.false. call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm, & urparm, fdum) write (iout,*) "Left SUMSL" penalty=v(10) iretcode=iv(1) nfun=iv(6) ngrad=iv(30) call targetfun(n,x,nf,f,uiparm,fdum) f0=f call scan(n,x,f,.false.) iscan=iscan+1 if (iv(1).eq.11 .and. cutoffviol) then write (iout,*) write (iout,'(a)') &"Revising the list of conformations because of cut-off violation." do iprot=1,nprot write (iout,'(a,f7.1)') & protname(iprot)(:ilen(protname(iprot))),enecut(iprot) do j=1,natlike(iprot)+2 write (iout,'(5(f7.1,f9.1,2x))') & (1.0d0/(Rgas*betaT(k,ibatch,iprot)), & Evar(j,k,iprot),k=1,nbeta(j,iprot)) enddo enddo #ifdef MPI cutoffeval=.false. call func1(n,5,x) #else call make_list(.true.,.false.,n,x) #endif goto 101 else if (iv(1).eq.11 .and. WhatsUp.eq.1) then write (iout,*) "Time up, terminating." goto 121 endif enddo else 102 continue cutoffviol=.false. cutoffeval=.false. write (iout,*) "MAXMIN",maxmin if (maxmin.gt.0) & call sumsl(n,d,x,targetfun,targetgrad,iv,liv,lv,v,uiparm, & urparm, fdum) write (iout,*) "Left SUMSL" penalty=v(10) iretcode=iv(1) nfun=nfun+iv(6) ngrad=ngrad+iv(30) niter=niter+iv(31) write (iout,*) "iretcode",iv(1)," cutoffviol",cutoffviol if (iv(1).eq.11 .and. cutoffviol) then write (iout,'(a)') &"Revising the list of conformations because of cut-off violation." do iprot=1,nprot write (iout,'(a,f7.1)') & protname(iprot)(:ilen(protname(iprot))),enecut(iprot) do j=1,natlike(iprot)+2 write (iout,'(5(f7.1,f9.1,2x))') & (1.0d0/(Rgas*betaT(k,ibatch,iprot)), & Evar(j,k,iprot),k=1,nbeta(j,iprot)) enddo enddo write (iout,'(a)') do iprot=1,nprot Evarmax(iprot)=Evar(1,1,iprot)*betaT(1,1,iprot) do ibatch=2,natlike(iprot)+2 do k=1,nbeta(ibatch,iprot) aux = Evar(k,ibatch,iprot)*betaT(k,ibatch,iprot) if (aux.gt.Evarmax(iprot)) Evarmax(iprot)=aux enddo enddo enddo do iprot=1,nprot aux=dmin1(0.5d0*enecut_min(iprot)+1.25d0*Evarmax(iprot) & ,enecut_max(iprot)) enecut(iprot)=dmax1(aux,enecut_min(iprot)) enddo write (iout,*) "Revised energy cutoffs" do iprot=1,nprot write (iout,*) "Protein",iprot,":",enecut(iprot) enddo call flush(iout) #ifdef MPI c c Send the order to the second master to revise the list of conformations. c cutoffeval=.false. call func1(n,5,x) #else call make_list(.true.,.false.,n,x) #endif call minimize_init(n,iv,v,d) iv(17)=iv(17)-nfun iv(18)=iv(18)-niter write (iout,*) "Diminishing IV(17) to",iv(17) write (iout,*) "Diminishing IV(18) to",iv(18) goto 102 else if (iv(1).eq.11 .and. WhatsUp.eq.1) then write (iout,*) "Time up, terminating." goto 121 endif endif 121 continue call targetfun(n,x,nf,f,uiparm,fdum) write (iout,*) write(iout,*) "Final function value:",f,v(10) write(iout,*) "Final value of the object function:",obf write(iout,*) "Number of target function evaluations:",nfun write(iout,*) "Number of target gradient evaluations:",ngrad write(iout,*) "SUMSL return code:",iretcode write(*,*) "Final function value:",f,v(10) write(*,*) "Final value of the object function:",obf write(*,*) "Number of target function evaluations:",nfun write(*,*) "Number of target gradient evaluations:",ngrad write(*,*) "SUMSL return code:",iretcode write (iout,*) call x2w(n,x) do i=1, n xm(i)=x(i) enddo open (88,file=restartname,status="unknown") do i=1,n write (88,*) x(i) enddo close(88) #ifdef CHECKGRAD write (iout,*) "Final checkgrad:" call checkgrad(n,x) #endif ELSE write (iout,*) "Unknown minimizer." ENDIF return end c----------------------------------------------------------------------------- real function funk(p) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" real p(max_paropt) double precision x(max_paropt) double precision f integer uiparm double precision obf double precision ufparm include "COMMON.IOUNITS" integer n,nf common /cc/ n integer i external fdum nf=nf+1 write (iout,*) "funk n",n," p",(p(i),i=1,n) do i=1,n x(i)=p(i) enddo c write (iout,*) "funk n",n," x",(x(i),i=1,n) call targetfun(n,x,nf,f,uiparm,fdum) funk=f return end c----------------------------------------------------------------------------- c logical function stopx(idum) c integer idum c stopx=.false. c return c end c----------------------------------------------------------------------------- subroutine checkgrad(n,x) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror integer status(MPI_STATUS_SIZE) #endif include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.OPTIM" include "COMMON.TIME1" include "COMMON.ENERGIES" include "COMMON.WEIGHTDER" include "COMMON.WEIGHTS" include "COMMON.CLASSES" include "COMMON.INTERACT" include "COMMON.NAMES" include "COMMON.VAR" include "COMMON.COMPAR" include "COMMON.TORSION" double precision energia(0:max_ene),epskl,eneps_num double precision f,fplus,fminus,fpot,fpotplus,xi,delta /5.0d-7/, & g(max_paropt),gnum(max_paropt),gfpot(max_paropt),gg(max_paropt), & gforce(max_paropt),rdum,zscder,pom,fchi,fchiplus, & sumlikder1(maxvar,maxT,maxprot) double precision chisquare_force integer n double precision x(n) integer uiparm,maxfun,i,ibatch,iprot,j,nf,k,kk,l,ii,iti,itj,ib, & l1,l2,jj,inat double precision nuave_orig(maxdimnat,maxT,maxnatlike,maxprot) #ifdef EPSCHECK double precision enetb_orig1(maxstr_proc,max_ene,maxprot) #endif integer icant external icant double precision rdif double precision urparm,aux,fac,eoldsc double precision & rder_num, & rder_all(max_paropt,maxdimnat,maxT,maxnatlike,maxprot), & zvtot_orig(maxT,maxprot),sumlik_orig(maxT,maxprot) external fdum integer ind_shield /25/ c write (iout,*) "Entered CHECKGRAD" call flush(iout) call targetfun(n,x,nf,f,uiparm,fdum) c write (iout,*) "f",f call targetgrad(n,x,nf,g,uiparm,rdum,fdum) c write (iout,*) "g",g(:n) write (iout,*) "Initial function & gradient evaluation completed." call flush(iout) #ifdef EPSCHECK if (mod_side) then do iprot=1,nprot do i=1,ntot(iprot) jj = i2ii(i,iprot) if (jj.gt.0) then do j=1,n_ene enetb_orig1(jj,j,iprot)=enetb(jj,j,iprot) enddo c write (iout,*) i,enetb_orig1(jj,1,iprot) endif enddo enddo write (iout,*) "Derivatives in epsilons" DO iprot=1,nprot call restore_molinfo(iprot) do i=1,ntot(iprot) jj = i2ii(i,iprot) if (jj.gt.0) then kk=0 c write (iout,*) "iprot",iprot," i",i call etotal(energia(0)) eoldsc=energia(1) do k=1,ntyp do l=1,k epskl=eps(k,l) eps(k,l)=eps(k,l)+delta eps(l,k)=eps(k,l) c write (iout,*) "***",k,l,epskl call prepare_sc kk=kk+1 call restore_ccoords(iprot,i) call etotal(energia(0)) eneps_num=(energia(1)-eoldsc) & /delta c write (iout,*) energia(1),enetb_orig1(jj,1,iprot) eps(k,l)=epskl eps(l,k)=epskl c write (iout,*) "---",k,l,epskl write(iout,'(2a4,3e20.10)') restyp(k),restyp(l), & eneps(kk,i,iprot),eneps_num, & (eneps_num-eneps(kk,i,iprot))*100 enddo enddo endif enddo ENDDO call prepare_sc endif #endif write (iout,*) "calling func1 0" call flush(iout) call func1(n,1,x) #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then fpot = rdif(n,x,gfpot,2) endif #endif c write (iout,*) "-----------------------------------" c write (iout,*) "Computing force chisquare" fchi = chisquare_force(n,x,gforce,2) c write (iout,*) "After first call to chisquare_force" c write (iout,*) "fchi",fchi c write (iout,*) "gforce",gforce(:n) c write (iout,*) "-----------------------------------" do iprot=1,nprot do inat=1,natlike(iprot) do ib=1,nbeta(inat+2,iprot) do i=1,natdim(inat,iprot) nuave_orig(i,ib,inat,iprot)=nuave(i,ib,inat,iprot) enddo enddo enddo enddo do iprot=1,nprot do ib=1,nbeta(1,iprot) sumlik_orig(ib,iprot)=sumlik(ib,iprot) enddo enddo do iprot=1,nprot do ib=1,nbeta(2,iprot) zvtot_orig(ib,iprot)=zvtot(ib,iprot) enddo enddo c write(iout,*) "exit func1 0" do iprot=1,nprot do ib=1,nbeta(2,iprot) do i=1,n zvdertot(i,ib,iprot)=0.0d0 enddo enddo enddo call zderiv ii=0 do k=1,n_ene if (imask(k).gt.0 .and. k.ne.ind_shield) ii=ii+1 enddo do iprot=1,nprot do ib=1,nbeta(1,iprot) do k=1,ii sumlikder1(k,ib,iprot)=sumlikder(k,ib,iprot) c write (iout,*) "ib",ib," k",k," sumlikder", c & sumlikder1(k,ib,iprot) enddo enddo enddo #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0 .and. torsion_pmf) then do i=1,nloctyp*(2*nloctyp-1) ii=ii+1 sumlikder1(ii,ib,iprot)=0.0d0 enddo endif #endif kk=ii do iprot=1,nprot c Maximum likelihood terms do ib=1,nbeta(1,iprot) kk=ii do i=1,nsingle_sc kk=kk+1 iti=ityp_ssc(i) sumlikder1(kk,ib,iprot)= & sumlikeps(icant(iti,iti),ib,iprot) do j=1,ntyp sumlikder1(kk,ib,iprot)= & sumlikder1(kk,ib,iprot)+ & sumlikeps(icant(iti,itj),ib,iprot) enddo enddo do i=1,npair_sc kk=kk+1 iti=ityp_psc(1,i) itj=ityp_psc(2,i) sumlikder1(kk,ib,iprot)= & sumlikeps(icant(iti,itj),ib,iprot) c write (iout,*) "kk",kk," iprot",iprot," ib",ib, c & " iti",iti," itj",itj, c & " ind",icant(iti,itj)," sumlikeps", c & sumlikeps(icant(iti,itj),ib,iprot) enddo do i=1,ntor_var kk=kk+1 sumlikder1(kk,ib,iprot)=sumliktor(i,ib,iprot) enddo do i=1,nang_var kk=kk+1 sumlikder1(kk,ib,iprot)=sumlikang(i,ib,iprot) enddo enddo c Heat capacity terms do ib=1,nbeta(2,iprot) kk=ii do i=1,nsingle_sc kk=kk+1 iti=ityp_ssc(i) zvdertot(kk,ib,iprot)= & zvepstot(icant(iti,iti),ib,iprot) do j=1,ntyp zvdertot(kk,ib,iprot)= & zvdertot(kk,ib,iprot)+ & zvepstot(icant(iti,j),ib,iprot) enddo enddo do i=1,npair_sc kk=kk+1 iti=ityp_psc(1,i) itj=ityp_psc(2,i) zvdertot(kk,ib,iprot)= & zvepstot(icant(iti,itj),ib,iprot) enddo do i=1,ntor_var kk=kk+1 zvdertot(kk,ib,iprot)=zvtortot(i,ib,iprot) enddo do i=1,nang_var kk=kk+1 zvdertot(kk,ib,iprot)=zvangtot(i,ib,iprot) enddo enddo enddo c Molecular properties do iprot=1,nprot do inat=1,natlike(iprot) do ib=1,nbeta(inat+2,iprot) call propderiv(ib,inat,iprot) do k=1,natdim(inat,iprot) do kk=1,ii rder_all(kk,k,ib,inat,iprot)=rder(kk,k) enddo kk=ii do i=1,nsingle_sc kk=kk+1 iti=ityp_ssc(i) rder_all(kk,k,ib,inat,iprot)= & reps(icant(iti,iti),k) do j=1,ntyp rder_all(kk,k,ib,inat,iprot)= & rder_all(kk,k,ib,inat,iprot)+ & reps(icant(iti,j),k) enddo enddo do i=1,npair_sc kk=kk+1 iti=ityp_psc(1,i) itj=ityp_psc(2,i) rder_all(kk,k,ib,inat,iprot)= & reps(icant(iti,itj),k) enddo kk=ii do i=1,ntor_var kk=kk+1 rder_all(kk,k,ib,inat,iprot)=rtor(i,k) enddo do i=1,nang_var kk=kk+1 rder_all(kk,k,ib,inat,iprot)=rang(i,k) enddo enddo enddo enddo enddo c Calculate numerical derivatives and compare with analytical derivatives do i=1,kk c write (iout,*) "gforce1:",gforce(:n) xi=x(i) x(i)=xi+xi*delta c write (iout,*) (x(k),k=1,kk) call func1(n,1,x) c write (iout,*) "after func1" c call flush(iout) #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then fpotplus = rdif(n,x,gfpot,1) c write (iout,*) "after rdif" c call flush(iout) endif #endif fchiplus = chisquare_force(n,x,gg,1) c write (iout,*) "gforce2:",gforce(:n) c write (iout,*) "Variable",i," fchiplus",fchiplus c Maximum likelihood do iprot=1,nprot write (iout,*) "Derivatives of maximum likelihood" do ib=1,nbeta(1,iprot) c write (iout,*) sumlik_orig(ib,iprot),sumlik(ib,iprot) pom=(sumlik(ib,iprot)-sumlik_orig(ib,iprot))/(xi*delta) write (iout,'(2i5,3e14.5)') iprot,ib, & sumlikder1(i,ib,iprot),pom,(pom-sumlikder1(i,ib,iprot)) & /(dabs(sumlikder1(i,ib,iprot))+1.0d-10) c write (iout,*) sumlik(ib,iprot), c & sumlik_orig(ib,iprot) enddo enddo c Heat capacity do iprot=1,nprot write (iout,*) "Derivatives of heat capacity" do ib=1,nbeta(2,iprot) pom=(zvtot(ib,iprot)-zvtot_orig(ib,iprot))/(xi*delta) write (iout,'(2i5,3e14.5)') iprot,ib, & zvdertot(i,ib,iprot),pom,(pom-zvdertot(i,ib,iprot)) & /(dabs(zvdertot(i,ib,iprot))+1.0d-10) c write (iout,*) zvtot(ib,ibatch,iprot), c & zvtot_orig(ib,ibatch,iprot) enddo enddo c Numerical derivatives of natlike properties write (iout,*) "Derivatives of molecular properties" do iprot=1,nprot do inat=1,natlike(iprot) do ib=1,nbeta(inat+2,iprot) do k=1,natdim(inat,iprot) rder_num=(nuave(k,ib,inat,iprot) & -nuave_orig(k,ib,inat,iprot))/(xi*delta) write (iout,'(4i3,2e14.5,f5.1)') iprot,inat,ib,k, & rder_all(kk,k,ib,inat,iprot), & rder_num,(rder_num-rder_all(kk,k,ib,inat,iprot))/ & (dabs(rder_all(kk,k,ib,inat,iprot))+1.0d-10)*100 enddo enddo enddo enddo #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then write(iout,*)"Derivative of fpot",(-fpotplus+fpot)/(xi*delta), & gfpot(i),(gfpot(i)-(-fpotplus+fpot)/(xi*delta))/gfpot(i) endif #endif c write (iout,*) "fchiplus",fchiplus," fchi",fchi," xi",xi, c & " delta",delta write (iout,*)"Derivatives of fchi",(-fchiplus+fchi)/(xi*delta), & gforce(i),(gforce(i)-(-fchiplus+fchi)/(xi*delta))/gforce(i) x(i)=xi enddo write (iout,*) do i=1,n xi=x(i) x(i)=xi+xi*delta c write (iout,*) "i",i," xplus ",x(i),xi call targetfun(n,x,nf,fplus,uiparm,fdum) x(i)=xi-xi*delta c write (iout,*) "i",i," xminus",x(i),xi call targetfun(n,x,nf,fminus,uiparm,fdum) x(i)=xi write (iout,*) "x",i,x(i)," fplus",fplus," fminus",fminus gnum(i)=0.5d0*(fplus-fminus)/(xi*delta) enddo write (iout,*) "Comparison of analytical and numerical gradient" do i=1,n write (iout,'(i5,2e15.7,f10.2)') i,g(i),gnum(i), & 1.0d2*(g(i)-gnum(i))/(dabs(g(i))+1.0d-10) enddo write (iout,*) return end c----------------------------------------------------------------------- subroutine prepare_sc implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.NAMES" include "COMMON.WEIGHTS" include "COMMON.INTERACT" include "COMMON.ENERGIES" include "COMMON.FFIELD" include "COMMON.TORSION" include "COMMON.IOUNITS" logical lprint integer i,j,ii,nvarr,iti,itj double precision rri double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2, & ratsig1,ratsig2,rsum_max,dwa16,r_augm lprint=.false. 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 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) 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) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa(i,j)=epsij*rrij*rrij bb(i,j)=-sigeps*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 return end