diff --git a/lib/reax/README b/lib/reax/README index 4345568624..c5dd5f0ae6 100644 --- a/lib/reax/README +++ b/lib/reax/README @@ -28,6 +28,7 @@ library source files and the LAMMPS pair_reax.cpp source file (from module REAX). It contains dimensions of statically-allocated arrays created by the ReaxFF library. The size of these arrays must be set small enough to avoid exceeding the available machine memory, and -large enough to fit the actual data generated by ReaxFF. If you -change the values in reax_defs.h, you must first rebuild the library -and then rebuild LAMMPS. +large enough to fit the actual data generated by ReaxFF. If you change +the values in reax_defs.h, you must first rebuild the library and then +rebuild LAMMPS. + diff --git a/lib/reax/cbka.blk b/lib/reax/cbka.blk index a5f3afb53a..4dbe0a36c6 100644 --- a/lib/reax/cbka.blk +++ b/lib/reax/cbka.blk @@ -54,6 +54,7 @@ $ elmol(nmolmax), $ elaf(nsort),vpq(nsort), $ rvdw(nsort),alf(nsort),eps(nsort),chat(nsort), + $ rcore(nsort,nsort),ecore(nsort,nsort),acore(nsort,nsort), $ vlp2(nsort), $ valp2(nsort),vincr(nsort), $ vval3(nsort), diff --git a/lib/reax/cbkc.blk b/lib/reax/cbkc.blk index a40dfb269f..2a2519eb2e 100644 --- a/lib/reax/cbkc.blk +++ b/lib/reax/cbkc.blk @@ -1,4 +1,5 @@ common - $/cbkc/ c(nat,3),cglobal(nattot,3),itag(nat) + $/cbkc/ c(nat,3),cglobal(nattot,3),itag(nat), + $chgglobal(nattot) diff --git a/lib/reax/reax_connect.F b/lib/reax/reax_connect.F index 09e48d3b5e..2b7ecc15c8 100644 --- a/lib/reax/reax_connect.F +++ b/lib/reax/reax_connect.F @@ -393,7 +393,6 @@ c$$$ call mdsav(1,qfile(nprob)) * lhb controls whether or not to unprune ghost bonds that * * may possibly form ghost hydrogen bonds. * * Setting it to 1 causes unpruning, and so is the safe option. * -* With the new midpoint method, lhb=0 may miss some h-bonds * * If lprune = 0, then pruning is not used, results are exact * * and lhb has no effect. * * * @@ -604,17 +603,19 @@ c but it does store bond entries in nubon2() ihhb2=nphb(ia(j2,1)) end if -* Only need to compute bonds where j3 is local +* Only need to compute bonds where j1 is local if (j1 .le. na_local) then if (ihhb1.eq.1.and.ihhb2.eq.2) then call dista2(j1,j2,dishb,dxm,dym,dzm) if (dishb.lt.hbcut) then do i23=1,ia(j1,2) !Search for acceptor atoms bound to donor atom - j3=ia(j1,2+i23) - if (nphb(ia(j3,1)).eq.2.and.j3.ne.j2) then - nboncol(nubon2(j1,i23))=-1 - ntmphb = ntmphb+1 + if (nboncol(nubon2(j1,i23)).eq.0) then + j3=ia(j1,2+i23) + if (nphb(ia(j3,1)).eq.2.and.j3.ne.j2) then + nboncol(nubon2(j1,i23))=-1 + ntmphb = ntmphb+1 + endif endif end do end if @@ -1475,7 +1476,7 @@ c ihhb2=nphb(ia(j2,1)) end if -* Only need to compute bonds where j3 is local +* Only need to compute bonds where j1 is local if (j1 .le. na_local) then if (ihhb1.eq.1.and.ihhb2.eq.2) then diff --git a/lib/reax/reax_inout.F b/lib/reax/reax_inout.F index 41cedba50a..f8675dce2a 100644 --- a/lib/reax/reax_inout.F +++ b/lib/reax/reax_inout.F @@ -33,6 +33,7 @@ #include "valang.blk" #include "cbksrtbon1.blk" #include "cbkchb.blk" + dimension rcore2(nsort),ecore2(nsort),acore2(nsort) ********************************************************************** * * * Read in force field * @@ -44,9 +45,6 @@ c$$$ write (65,*) 'In ffinpt' c$$$ call timer(65) c$$$ close (65) c$$$ end if - - open (4,file='ffield.reax',status='old') - rewind (4) iline=0 read (4,'(a40)',end=990,err=990)qffield @@ -96,7 +94,7 @@ c$$$ end if $bo131(i1),bo132(i1),bo133(i1),sigqeq(i1),default iline=iline+1 read (4,1250,end=990,err=990)vovun(i1),vval1(i1),vrom, - $vval3(i1),vval4(i1) + $vval3(i1),vval4(i1),rcore2(i1),ecore2(i1),acore2(i1) iline=iline+1 idef(i1)=int(default) nphb(i1)=int(vnphb) @@ -108,6 +106,9 @@ c$$$ end if ********************************************************************** do i1=1,nso do i2=1,nso + rcore(i1,i2)=sqrt(rcore2(i1)*rcore2(i2)) + ecore(i1,i2)=sqrt(ecore2(i1)*ecore2(i2)) + acore(i1,i2)=sqrt(acore2(i1)*acore2(i2)) p1co(i1,i2)=sqrt(4.0*rvdw(i1)*rvdw(i2)) p2co(i1,i2)=sqrt(eps(i1)*eps(i2)) p3co(i1,i2)=sqrt(alf(i1)*alf(i2)) @@ -261,6 +262,7 @@ c$$$ end if 990 write (*,*)'Error or end-of-file reading unit 4 on line:',iline stop 999 continue + close(4) ********************************************************************** * * * Format part * @@ -2182,7 +2184,8 @@ c$$$ end if $,end=40,err=40) $ir,qlabel(na+1),qresi1(na+1),qresi2(na+1),qresi3(na+1), $cglobal(na+1,1),cglobal(na+1,2), - $cglobal(na+1,3),qffty(na+1),ibgr1(na+1),ibgr2(na+1),chgbgf(na+1) + $cglobal(na+1,3),qffty(na+1),ibgr1(na+1),ibgr2(na+1), + $chgglobal(na+1) else stop 'Unsupported Biograf-version' end if diff --git a/lib/reax/reax_poten.F b/lib/reax/reax_poten.F index 62e810d47f..a407a0b8ca 100644 --- a/lib/reax/reax_poten.F +++ b/lib/reax/reax_poten.F @@ -2805,8 +2805,6 @@ c$$$ end if exphuc=exp(-vpar(24)*boc) bocor4=(1.0-exphua)*(1.0-exphub)*(1.0-exphuc) eth=hsin*ethhulp*bocor4 - estrain(j(2))=estrain(j(2))+0.50*eth - estrain(j(3))=estrain(j(3))+0.50*eth detdar=hsin*bocor4*(0.50*v1(ity)-2.0*v2(ity)*bocor2*arg+ $v3(ity)*(6.0*arg2-1.5d0)) diff --git a/lib/reax/reax_reac.F b/lib/reax/reax_reac.F index c638944649..cdac3e2546 100644 --- a/lib/reax/reax_reac.F +++ b/lib/reax/reax_reac.F @@ -93,9 +93,10 @@ c$$$ end if call lonpar call covbon call ovcor + call srtang !Determine valency angles call srttor !Determine torsion angles - call srtoop !Determine out of plane angles +* call srtoop !Determine out of plane angles call srthb !Determine hydrogen bonds call calval