diff --git a/source/analysis.f b/source/analysis.f index 5db8e88c..18e761cc 100644 --- a/source/analysis.f +++ b/source/analysis.f @@ -26,6 +26,7 @@ subroutine analysis (energy) use iounit use limits use potent + use reppot use vdwpot implicit none integer i @@ -227,7 +228,11 @@ subroutine analysis (energy) if (vdwtyp .eq. 'BUFFERED-14-7') call ehal3 if (vdwtyp .eq. 'GAUSSIAN') call egauss3 end if - if (use_repel) call erepel3 + if (use_repel) then + call erepel3 + else if (use_xrepel) then + call exrepel3 + end if if (use_disp) call edisp3 c c call any miscellaneous energy component routines diff --git a/source/analyze.f b/source/analyze.f index 4cfb9f7b..8167a2cf 100644 --- a/source/analyze.f +++ b/source/analyze.f @@ -384,7 +384,7 @@ subroutine systyze c c details of vdw potential energy functional form c - if (use_vdw .or. use_repel .or. use_disp) then + if (use_vdw .or. use_repel .or. use_xrepel .or. use_disp) then write (iout,70) 70 format () end if @@ -413,20 +413,25 @@ subroutine systyze call justify (value) write (iout,100) value 100 format (' VDW Function',16x,a20) + else if (use_xrepel) then + value = 'EXCHANGE REPULSION' + call justify (value) + write (iout,110) value + 110 format (' VDW Function',16x,a20) end if if (use_disp) then value = 'DAMPED DISPERSION' call justify (value) - write (iout,110) value - 110 format (' VDW Function',16x,a20) + write (iout,120) value + 120 format (' VDW Function',16x,a20) end if c c details of dispersion particle mesh Ewald calculation c if (use_dewald) then - write (iout,120) adewald,dewaldcut,ndfft1, + write (iout,130) adewald,dewaldcut,ndfft1, & ndfft2,ndfft3,bsdorder - 120 format (/,' PME for Dispersion :', + 130 format (/,' PME for Dispersion :', & //,' Ewald Coefficient',19x,f12.4, & /,' Real-Space Cutoff',19x,f12.4, & /,' Grid Dimensions',21x,3i4, @@ -436,59 +441,59 @@ subroutine systyze c details of electrostatic energy functional form c if (use_charge .or. use_dipole .or. use_mpole .or. use_polar) then - write (iout,130) - 130 format () + write (iout,140) + 140 format () end if if (use_charge) then value = 'PARTIAL CHARGE' call justify (value) - write (iout,140) value - 140 format (' Electrostatics',14x,a20) + write (iout,150) value + 150 format (' Electrostatics',14x,a20) end if if (use_dipole) then value = 'BOND DIPOLE' call justify (value) - write (iout,150) value - 150 format (' Electrostatics',14x,a20) + write (iout,160) value + 160 format (' Electrostatics',14x,a20) end if if (use_mpole) then value = 'ATOMIC MULTIPOLE' call justify (value) - write (iout,160) value - 160 format (' Electrostatics',14x,a20) + write (iout,170) value + 170 format (' Electrostatics',14x,a20) end if if (use_chgpen) then value = 'CHARGE PENETRATION' call justify (value) - write (iout,170) value - 170 format (' Electrostatics',14x,a20) + write (iout,180) value + 180 format (' Electrostatics',14x,a20) end if if (use_chgtrn) then value = 'CHARGE TRANSFER' call justify (value) - write (iout,180) value - 180 format (' Induction',19x,a20) + write (iout,190) value + 190 format (' Induction',19x,a20) end if if (use_polar) then value = 'INDUCED DIPOLE' call justify (value) - write (iout,190) value - 190 format (' Induction',19x,a20) + write (iout,200) value + 200 format (' Induction',19x,a20) value = poltyp call justify (value) - write (iout,200) value - 200 format (' Polarization',16x,a20) + write (iout,210) value + 210 format (' Polarization',16x,a20) if (use_thole) then value = 'THOLE DAMPING' call justify (value) - write (iout,210) value - 210 format (' Polarization',16x,a20) + write (iout,220) value + 220 format (' Polarization',16x,a20) end if if (use_chgpen) then value = 'CHGPEN DAMPING' call justify (value) - write (iout,220) value - 220 format (' Polarization',16x,a20) + write (iout,230) value + 230 format (' Polarization',16x,a20) end if end if c @@ -497,9 +502,9 @@ subroutine systyze if (use_ewald) then value = boundary call justify (value) - write (iout,230) aeewald,ewaldcut,nefft1,nefft2, + write (iout,240) aeewald,ewaldcut,nefft1,nefft2, & nefft3,bseorder,value - 230 format (/,' PME for Electrostatics :', + 240 format (/,' PME for Electrostatics :', & //,' Ewald Coefficient',19x,f12.4, & /,' Real-Space Cutoff',19x,f12.4, & /,' Grid Dimensions',21x,3i4, @@ -663,6 +668,13 @@ subroutine partyze fstr = '('' Repulsion'',18x,'//form2//')' end if write (iout,fstr) er,ner + else if (use_xrepel .and. (ner.ne.0.or.er.ne.0.0d0)) then + if (abs(er) .lt. 1.0d10) then + fstr = '('' Repulsion'',18x,'//form1//')' + else + fstr = '('' Repulsion'',18x,'//form2//')' + end if + write (iout,fstr) er,ner end if if (use_disp .and. (nedsp.ne.0.or.edsp.ne.0.0d0)) then fstr = '('' Dispersion'',17x,'//form1//')' @@ -1226,6 +1238,7 @@ subroutine paramyze (active) use urey use vdw use vdwpot + use xrepel implicit none integer i,j,k integer ia,ib,ic @@ -1314,50 +1327,53 @@ subroutine paramyze (active) if (use_repel .and. nrep.ne.0) then write (iout,180) nrep 180 format (' Repulsion Sites',18x,i15) + else if (use_xrepel .and. nxrep.ne.0) then + write (iout,190) nxrep + 190 format (' Repulsion Sites',18x,i15) end if if (use_disp .and. ndisp.ne.0) then - write (iout,190) ndisp - 190 format (' Dispersion Sites',17x,i15) + write (iout,200) ndisp + 200 format (' Dispersion Sites',17x,i15) end if if (use_charge .and. nion.ne.0) then - write (iout,200) nion - 200 format (' Atomic Partial Charges',11x,i15) + write (iout,210) nion + 210 format (' Atomic Partial Charges',11x,i15) end if if (use_dipole .and. ndipole.ne.0) then - write (iout,210) ndipole - 210 format (' Bond Dipole Moments',14x,i15) + write (iout,220) ndipole + 220 format (' Bond Dipole Moments',14x,i15) end if if (use_mpole .and. npole.ne.0) then - write (iout,220) npole - 220 format (' Atomic Multipoles',16x,i15) + write (iout,230) npole + 230 format (' Atomic Multipoles',16x,i15) end if if (use_chgpen .and. ncp.ne.0) then - write (iout,230) ncp - 230 format (' Charge Penetration',15x,i15) + write (iout,240) ncp + 240 format (' Charge Penetration',15x,i15) end if if (use_polar .and. npolar.ne.0) then - write (iout,240) npolar - 240 format (' Polarizable Sites',16x,i15) + write (iout,250) npolar + 250 format (' Polarizable Sites',16x,i15) end if if (use_chgtrn .and. nct.ne.0) then - write (iout,250) nct - 250 format (' Charge Transfer Sites',12x,i15) + write (iout,260) nct + 260 format (' Charge Transfer Sites',12x,i15) end if if (use_chgflx .and. nbflx.ne.0) then - write (iout,260) nbflx - 260 format (' Bond Charge Flux',17x,i15) + write (iout,270) nbflx + 270 format (' Bond Charge Flux',17x,i15) end if if (use_chgflx .and. naflx.ne.0) then - write (iout,270) naflx - 270 format (' Angle Charge Flux',16x,i15) + write (iout,280) naflx + 280 format (' Angle Charge Flux',16x,i15) end if if (use_orbit .and. norbit.ne.0) then - write (iout,280) norbit - 280 format (' Pisystem Atoms',19x,i15) + write (iout,290) norbit + 290 format (' Pisystem Atoms',19x,i15) end if if (use_orbit .and. nbpi.ne.0) then - write (iout,290) nbpi - 290 format (' Conjugated Pi-Bonds',14x,i15) + write (iout,300) nbpi + 300 format (' Conjugated Pi-Bonds',14x,i15) end if c c parameters used for molecular mechanics atom types @@ -1367,15 +1383,15 @@ subroutine paramyze (active) if (active(i)) then if (header) then header = .false. - write (iout,300) - 300 format (/,' Atom Definition Parameters :', + write (iout,310) + 310 format (/,' Atom Definition Parameters :', & //,3x,'Atom',2x,'Symbol',2x,'Type', & 2x,'Class',2x,'Atomic',3x,'Mass', & 2x,'Valence',2x,'Description',/) end if - write (iout,310) i,name(i),type(i),class(i),atomic(i), + write (iout,320) i,name(i),type(i),class(i),atomic(i), & mass(i),valnum(i),story(i) - 310 format (i6,5x,a3,2i7,i6,f10.3,i5,5x,a24) + 320 format (i6,5x,a3,2i7,i6,f10.3,i5,5x,a24) end if end do c @@ -1389,12 +1405,12 @@ subroutine paramyze (active) if (active(ia) .or. active(ib)) then if (header) then header = .false. - write (iout,320) - 320 format (/,' Bond Stretching Parameters :', + write (iout,330) + 330 format (/,' Bond Stretching Parameters :', & //,10x,'Atom Numbers',25x,'KS',7x,'Bond',/) end if - write (iout,330) i,ia,ib,bk(i),bl(i) - 330 format (i6,3x,2i6,19x,f10.3,f10.4) + write (iout,340) i,ia,ib,bk(i),bl(i) + 340 format (i6,3x,2i6,19x,f10.3,f10.4) end if end do end if @@ -1410,23 +1426,23 @@ subroutine paramyze (active) if (active(ia) .or. active(ib) .or. active(ic)) then if (header) then header = .false. - write (iout,340) - 340 format (/,' Angle Bending Parameters :', + write (iout,350) + 350 format (/,' Angle Bending Parameters :', & //,13x,'Atom Numbers',22x,'KB', & 6x,'Angle',3x,'Fold',4x,'Type',/) end if if (angtyp(i) .eq. 'HARMONIC') then - write (iout,350) i,ia,ib,ic,ak(i),anat(i) - 350 format (i6,3x,3i6,13x,2f10.3) - else if (angtyp(i) .eq. 'IN-PLANE') then write (iout,360) i,ia,ib,ic,ak(i),anat(i) - 360 format (i6,3x,3i6,13x,2f10.3,9x,'In-Plane') - else if (angtyp(i) .eq. 'LINEAR') then + 360 format (i6,3x,3i6,13x,2f10.3) + else if (angtyp(i) .eq. 'IN-PLANE') then write (iout,370) i,ia,ib,ic,ak(i),anat(i) - 370 format (i6,3x,3i6,13x,2f10.3,9x,'Linear') + 370 format (i6,3x,3i6,13x,2f10.3,9x,'In-Plane') + else if (angtyp(i) .eq. 'LINEAR') then + write (iout,380) i,ia,ib,ic,ak(i),anat(i) + 380 format (i6,3x,3i6,13x,2f10.3,9x,'Linear') else if (angtyp(i) .eq. 'FOURIER ') then - write (iout,380) i,ia,ib,ic,ak(i),anat(i),afld(i) - 380 format (i6,3x,3i6,13x,2f10.3,f7.1,2x,'Fourier') + write (iout,390) i,ia,ib,ic,ak(i),anat(i),afld(i) + 390 format (i6,3x,3i6,13x,2f10.3,f7.1,2x,'Fourier') end if end if end do @@ -1444,8 +1460,8 @@ subroutine paramyze (active) if (active(ia) .or. active(ib) .or. active(ic)) then if (header) then header = .false. - write (iout,390) - 390 format (/,' Stretch-Bend Parameters :', + write (iout,400) + 400 format (/,' Stretch-Bend Parameters :', & //,13x,'Atom Numbers',8x,'KSB 1',5x,'KSB 2', & 6x,'Angle',3x,'Bond 1',3x,'Bond 2',/) end if @@ -1453,9 +1469,9 @@ subroutine paramyze (active) blc = 0.0d0 if (isb(2,i) .ne. 0) bla = bl(isb(2,i)) if (isb(3,i) .ne. 0) blc = bl(isb(3,i)) - write (iout,400) i,ia,ib,ic,sbk(1,i),sbk(2,i), + write (iout,410) i,ia,ib,ic,sbk(1,i),sbk(2,i), & anat(k),bla,blc - 400 format (i6,3x,3i6,1x,2f10.3,2x,f9.3,2f9.4) + 410 format (i6,3x,3i6,1x,2f10.3,2x,f9.3,2f9.4) end if end do end if @@ -1471,13 +1487,13 @@ subroutine paramyze (active) if (active(ia) .or. active(ic)) then if (header) then header = .false. - write (iout,410) - 410 format (/,' Urey-Bradley Parameters :', + write (iout,420) + 420 format (/,' Urey-Bradley Parameters :', & //,13x,'Atom Numbers',21x,'KUB', & 4x,'Distance',/) end if - write (iout,420) i,ia,ib,ic,uk(i),ul(i) - 420 format (i6,3x,3i6,13x,f10.3,f10.4) + write (iout,430) i,ia,ib,ic,uk(i),ul(i) + 430 format (i6,3x,3i6,13x,f10.3,f10.4) end if end do end if @@ -1496,12 +1512,12 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,430) - 430 format (/,' Out-of-Plane Bend Parameters :', + write (iout,440) + 440 format (/,' Out-of-Plane Bend Parameters :', & //,17x,'Atom Numbers',19x,'KOPB',/) end if - write (iout,440) i,id,ib,ia,ic,opbk(i) - 440 format (i6,3x,4i6,9x,f10.3) + write (iout,450) i,id,ib,ia,ic,opbk(i) + 450 format (i6,3x,4i6,9x,f10.3) end if end do end if @@ -1519,12 +1535,12 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,450) - 450 format (/,' Out-of-Plane Distance Parameters :', + write (iout,460) + 460 format (/,' Out-of-Plane Distance Parameters :', & //,17x,'Atom Numbers',19x,'KOPD',/) end if - write (iout,460) i,ia,ib,ic,id,opdk(i) - 460 format (i6,3x,4i6,9x,f10.3) + write (iout,470) i,ia,ib,ic,id,opdk(i) + 470 format (i6,3x,4i6,9x,f10.3) end if end do end if @@ -1542,13 +1558,13 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,470) - 470 format (/,' Improper Dihedral Parameters :', + write (iout,480) + 480 format (/,' Improper Dihedral Parameters :', & //,17x,'Atom Numbers',19x,'KID', & 4x,'Dihedral',/) end if - write (iout,480) i,ia,ib,ic,id,kprop(i),vprop(i) - 480 format (i6,3x,4i6,9x,2f10.4) + write (iout,490) i,ia,ib,ic,id,kprop(i),vprop(i) + 490 format (i6,3x,4i6,9x,2f10.4) end if end do end if @@ -1566,8 +1582,8 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,490) - 490 format (/,' Improper Torsion Parameters :', + write (iout,500) + 500 format (/,' Improper Torsion Parameters :', & //,17x,'Atom Numbers',11x, & 'Amplitude, Phase and Periodicity',/) end if @@ -1591,20 +1607,20 @@ subroutine paramyze (active) phase(j) = itors3(2,i) end if if (j .eq. 0) then - write (iout,500) i,ia,ib,ic,id - 500 format (i6,3x,4i6) + write (iout,510) i,ia,ib,ic,id + 510 format (i6,3x,4i6) else if (j .eq. 1) then - write (iout,510) i,ia,ib,ic,id, + write (iout,520) i,ia,ib,ic,id, & ampli(1),phase(1),fold(1) - 510 format (i6,3x,4i6,10x,f10.3,f8.1,i4) + 520 format (i6,3x,4i6,10x,f10.3,f8.1,i4) else if (j .eq. 2) then - write (iout,520) i,ia,ib,ic,id,(ampli(k), + write (iout,530) i,ia,ib,ic,id,(ampli(k), & phase(k),fold(k),k=1,j) - 520 format (i6,3x,4i6,2x,2(f10.3,f6.1,i4)) + 530 format (i6,3x,4i6,2x,2(f10.3,f6.1,i4)) else - write (iout,530) i,ia,ib,ic,id,(ampli(k), + write (iout,540) i,ia,ib,ic,id,(ampli(k), & nint(phase(k)),fold(k),k=1,j) - 530 format (i6,3x,4i6,4x,3(f8.3,i4,'/',i1)) + 540 format (i6,3x,4i6,4x,3(f8.3,i4,'/',i1)) end if end if end do @@ -1623,8 +1639,8 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,540) - 540 format (/,' Torsional Angle Parameters :', + write (iout,550) + 550 format (/,' Torsional Angle Parameters :', & //,17x,'Atom Numbers',11x, & 'Amplitude, Phase and Periodicity',/) end if @@ -1666,12 +1682,12 @@ subroutine paramyze (active) phase(j) = tors6(2,i) end if if (j .eq. 0) then - write (iout,550) i,ia,ib,ic,id - 550 format (i6,3x,4i6) + write (iout,560) i,ia,ib,ic,id + 560 format (i6,3x,4i6) else - write (iout,560) i,ia,ib,ic,id,(ampli(k), + write (iout,570) i,ia,ib,ic,id,(ampli(k), & nint(phase(k)),fold(k),k=1,j) - 560 format (i6,3x,4i6,4x,6(f8.3,i4,'/',i1)) + 570 format (i6,3x,4i6,4x,6(f8.3,i4,'/',i1)) end if end if end do @@ -1692,12 +1708,12 @@ subroutine paramyze (active) & active(id) .or. active(ie) .or. active(ig)) then if (header) then header = .false. - write (iout,570) - 570 format (/,' Pi-Orbital Torsion Parameters :', + write (iout,580) + 580 format (/,' Pi-Orbital Torsion Parameters :', & //,10x,'Atom Numbers',19x,'Amplitude',/) end if - write (iout,580) i,ic,id,kpit(i) - 580 format (i6,3x,2i6,19x,f10.4) + write (iout,590) i,ic,id,kpit(i) + 590 format (i6,3x,2i6,19x,f10.4) end if end do end if @@ -1716,8 +1732,8 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,590) - 590 format (/,' Stretch-Torsion Parameters :', + write (iout,600) + 600 format (/,' Stretch-Torsion Parameters :', & //,17x,'Atom Numbers',10x,'Bond', & 5x,'Amplitude and Phase (1-3 Fold)',/) end if @@ -1739,11 +1755,11 @@ subroutine paramyze (active) phase(8) = tors2(2,k) ampli(9) = kst(9,i) phase(9) = tors3(2,k) - write (iout,600) i,ia,ib,ic,id, + write (iout,610) i,ia,ib,ic,id, & '1st',(ampli(k),nint(phase(k)),k=1,3), & '2nd',(ampli(k),nint(phase(k)),k=4,6), & '3rd',(ampli(k),nint(phase(k)),k=7,9) - 600 format (i6,3x,4i6,7x,a3,3x,3(f7.3,i4), + 610 format (i6,3x,4i6,7x,a3,3x,3(f7.3,i4), & /,40x,a3,3x,3(f7.3,i4), & /,40x,a3,3x,3(f7.3,i4)) end if @@ -1764,8 +1780,8 @@ subroutine paramyze (active) & active(ic) .or. active(id)) then if (header) then header = .false. - write (iout,610) - 610 format (/,' Angle-Torsion Parameters :', + write (iout,620) + 620 format (/,' Angle-Torsion Parameters :', & //,17x,'Atom Numbers',10x,'Angle', & 4x,'Amplitude and Phase (1-3 Fold)',/) end if @@ -1781,10 +1797,10 @@ subroutine paramyze (active) phase(5) = tors2(2,k) ampli(6) = kant(6,i) phase(6) = tors3(2,k) - write (iout,620) i,ia,ib,ic,id, + write (iout,630) i,ia,ib,ic,id, & '1st',(ampli(k),nint(phase(k)),k=1,3), & '2nd',(ampli(k),nint(phase(k)),k=4,6) - 620 format (i6,3x,4i6,7x,a3,3x,3(f7.3,i4), + 630 format (i6,3x,4i6,7x,a3,3x,3(f7.3,i4), & /,40x,a3,3x,3(f7.3,i4)) end if end do @@ -1805,14 +1821,14 @@ subroutine paramyze (active) & active(id) .or. active(ie)) then if (header) then header = .false. - write (iout,630) - 630 format (/,' Torsion-Torsion Parameters :', + write (iout,640) + 640 format (/,' Torsion-Torsion Parameters :', & //,20x,'Atom Numbers',15x,'Spline Grid', & 6x,'Tier',/) end if j = itt(2,i) - write (iout,640) i,ia,ib,ic,id,ie,tnx(j),tny(j),ttier(j) - 640 format (i6,3x,5i6,7x,2i6,7x,a3) + write (iout,650) i,ia,ib,ic,id,ie,tnx(j),tny(j),ttier(j) + 650 format (i6,3x,5i6,7x,2i6,7x,a3) end if end do end if @@ -1827,8 +1843,8 @@ subroutine paramyze (active) k = k + 1 if (header) then header = .false. - write (iout,650) - 650 format (/,' Van der Waals Parameters :', + write (iout,660) + 660 format (/,' Van der Waals Parameters :', & //,10x,'Atom Number',7x,'Size', & 3x,'Epsilon',3x,'Size 1-4', & 3x,'Eps 1-4',3x,'Reduction',/) @@ -1840,11 +1856,11 @@ subroutine paramyze (active) if (radsiz .eq. 'DIAMETER') radj = 2.0d0 * radj if (radtyp .eq. 'SIGMA') radj = radj / twosix if (reduct(j) .eq. 0.0d0) then - write (iout,660) k,i,radj,eps(j) - 660 format (i6,3x,i6,7x,2f10.4) + write (iout,670) k,i,radj,eps(j) + 670 format (i6,3x,i6,7x,2f10.4) else - write (iout,670) k,i,radj,eps(j),reduct(j) - 670 format (i6,3x,i6,7x,2f10.4,22x,f10.4) + write (iout,680) k,i,radj,eps(j),reduct(j) + 680 format (i6,3x,i6,7x,2f10.4,22x,f10.4) end if else radj = rad(j) @@ -1858,12 +1874,12 @@ subroutine paramyze (active) rad4j = rad4j / twosix end if if (reduct(j) .eq. 0.0d0) then - write (iout,680) k,i,radj,eps(j),rad4j,eps4(j) - 680 format (i6,3x,i6,7x,2f10.4,1x,2f10.4) + write (iout,690) k,i,radj,eps(j),rad4j,eps4(j) + 690 format (i6,3x,i6,7x,2f10.4,1x,2f10.4) else - write (iout,690) k,i,radj,eps(j),rad4j, + write (iout,700) k,i,radj,eps(j),rad4j, & eps4(j),reduct(j) - 690 format (i6,3x,i6,7x,2f10.4,1x,2f10.4,1x,f10.4) + 700 format (i6,3x,i6,7x,2f10.4,1x,2f10.4,1x,f10.4) end if end if end if @@ -1879,13 +1895,32 @@ subroutine paramyze (active) if (active(ia)) then if (header) then header = .false. - write (iout,700) - 700 format (/,' Pauli Repulsion Parameters :', + write (iout,710) + 710 format (/,' Pauli Repulsion Parameters :', & //,10x,'Atom Number',25x,'Size',6x,'Damp', & 3x,'Valence',/) end if - write (iout,710) i,ia,sizpr(i),dmppr(i),elepr(i) - 710 format (i6,3x,i6,25x,2f10.4,f10.3) + write (iout,720) i,ia,sizpr(i),dmppr(i),elepr(i) + 720 format (i6,3x,i6,25x,2f10.4,f10.3) + end if + end do +c +c parameters used for exchange repulsion interactions +c + else if (use_xrepel) then + header = .true. + do i = 1, nxrep + ia = ixrep(i) + if (active(ia)) then + if (header) then + header = .false. + write (iout,730) + 730 format (/,' Exchange Repulsion Parameters :', + & //,10x,'Atom Number',23x,'Charge',6x,'Damp', + & 2x,'PS-Ratio',/) + end if + write (iout,740) i,ia,zpxr(i),dmppxr(i),crpxr(i) + 740 format (i6,3x,i6,25x,2f10.4,f10.3) end if end do end if @@ -1899,12 +1934,12 @@ subroutine paramyze (active) if (active(ia)) then if (header) then header = .false. - write (iout,720) - 720 format (/,' Damped Dispersion Parameters :', + write (iout,750) + 750 format (/,' Damped Dispersion Parameters :', & //,10x,'Atom Number',26x,'C6',7x,'Damp',/) end if - write (iout,730) i,ia,csix(i),adisp(i) - 730 format (i6,3x,i6,25x,f10.3,f10.4) + write (iout,760) i,ia,csix(i),adisp(i) + 760 format (i6,3x,i6,25x,f10.3,f10.4) end if end do end if @@ -1920,18 +1955,18 @@ subroutine paramyze (active) if (active(ia) .or. active(ic)) then if (header) then header = .false. - write (iout,740) - 740 format (/,' Atomic Partial Charge Parameters :', + write (iout,770) + 770 format (/,' Atomic Partial Charge Parameters :', & /,45x,'Neighbor',3x,'Cutoff', & /,10x,'Atom Number',13x,'Charge', & 7x,'Site',6x,'Site',/) end if if (ia.eq.ib .and. ia.eq.ic) then - write (iout,750) i,ia,pchg(ia) - 750 format (i6,3x,i6,15x,f10.4) + write (iout,780) i,ia,pchg(ia) + 780 format (i6,3x,i6,15x,f10.4) else - write (iout,760) i,ia,pchg(ia),ib,ic - 760 format (i6,3x,i6,15x,f10.4,5x,i6,4x,i6) + write (iout,790) i,ia,pchg(ia),ib,ic + 790 format (i6,3x,i6,15x,f10.4,5x,i6,4x,i6) end if end if end do @@ -1947,13 +1982,13 @@ subroutine paramyze (active) if (active(ia) .or. active(ib)) then if (header) then header = .false. - write (iout,770) - 770 format (/,' Bond Dipole Moment Parameters :', + write (iout,800) + 800 format (/,' Bond Dipole Moment Parameters :', & //,10x,'Atom Numbers',22x,'Dipole', & 3x,'Position',/) end if - write (iout,780) i,ia,ib,bdpl(i),sdpl(i) - 780 format (i6,3x,2i6,19x,f10.4,f10.3) + write (iout,810) i,ia,ib,bdpl(i),sdpl(i) + 810 format (i6,3x,2i6,19x,f10.4,f10.3) end if end do end if @@ -1967,8 +2002,8 @@ subroutine paramyze (active) if (active(ia)) then if (header) then header = .false. - write (iout,790) - 790 format (/,' Atomic Multipole Parameters :', + write (iout,820) + 820 format (/,' Atomic Multipole Parameters :', & //,11x,'Atom',3x,'Z-Axis',1x,'X-Axis', & 1x,'Y-Axis',2x, & 'Frame',11x,'Multipole Moments',/) @@ -1985,28 +2020,28 @@ subroutine paramyze (active) mpl(j) = 3.0d0 * pole(j,ia) / bohr**2 end do if (izaxe .eq. 0) then - write (iout,800) i,ia,polaxe(i), + write (iout,830) i,ia,polaxe(i), & (mpl(j),j=1,5),mpl(8),mpl(9), & (mpl(j),j=11,13) - 800 format (i6,3x,i6,25x,a8,2x,f9.5,/,50x,3f9.5, + 830 format (i6,3x,i6,25x,a8,2x,f9.5,/,50x,3f9.5, & /,50x,f9.5,/,50x,2f9.5,/,50x,3f9.5) else if (ixaxe .eq. 0) then - write (iout,810) i,ia,izaxe,polaxe(i), + write (iout,840) i,ia,izaxe,polaxe(i), & (mpl(j),j=1,5),mpl(8),mpl(9), & (mpl(j),j=11,13) - 810 format (i6,3x,i6,1x,i7,17x,a8,2x,f9.5,/,50x,3f9.5, + 840 format (i6,3x,i6,1x,i7,17x,a8,2x,f9.5,/,50x,3f9.5, & /,50x,f9.5,/,50x,2f9.5,/,50x,3f9.5) else if (iyaxe .eq. 0) then - write (iout,820) i,ia,izaxe,ixaxe,polaxe(i), + write (iout,850) i,ia,izaxe,ixaxe,polaxe(i), & (mpl(j),j=1,5),mpl(8),mpl(9), & (mpl(j),j=11,13) - 820 format (i6,3x,i6,1x,2i7,10x,a8,2x,f9.5,/,50x,3f9.5, + 850 format (i6,3x,i6,1x,2i7,10x,a8,2x,f9.5,/,50x,3f9.5, & /,50x,f9.5,/,50x,2f9.5,/,50x,3f9.5) else - write (iout,830) i,ia,izaxe,ixaxe,iyaxe,polaxe(i), + write (iout,860) i,ia,izaxe,ixaxe,iyaxe,polaxe(i), & (mpl(j),j=1,5),mpl(8),mpl(9), & (mpl(j),j=11,13) - 830 format (i6,3x,i6,1x,3i7,3x,a8,2x,f9.5,/,50x,3f9.5, + 860 format (i6,3x,i6,1x,3i7,3x,a8,2x,f9.5,/,50x,3f9.5, & /,50x,f9.5,/,50x,2f9.5,/,50x,3f9.5) end if end if @@ -2022,13 +2057,13 @@ subroutine paramyze (active) if (active(ia)) then if (header) then header = .false. - write (iout,840) - 840 format (/,' Charge Penetration Parameters :', + write (iout,870) + 870 format (/,' Charge Penetration Parameters :', & //,10x,'Atom Number',25x,'Core',3x,'Valence', & 6x,'Damp',/) end if - write (iout,850) i,ia,pcore(ia),pval(ia),palpha(ia) - 850 format (i6,3x,i6,25x,2f10.3,f10.4) + write (iout,880) i,ia,pcore(ia),pval(ia),palpha(ia) + 880 format (i6,3x,i6,25x,2f10.3,f10.4) end if end do end if @@ -2043,34 +2078,34 @@ subroutine paramyze (active) if (header) then header = .false. if (use_tholed) then - write (iout,860) - 860 format (/,' Dipole Polarizability Parameters :', + write (iout,890) + 890 format (/,' Dipole Polarizability Parameters :', & //,10x,'Atom Number',5x,'Alpha',4x,'Thole', & 3x,'TholeD',6x,'Polarization Group',/) else if (use_thole) then - write (iout,870) - 870 format (/,' Dipole Polarizability Parameters :', + write (iout,900) + 900 format (/,' Dipole Polarizability Parameters :', & //,10x,'Atom Number',5x,'Alpha',4x,'Thole', & 6x,'Polarization Group',/) else - write (iout,880) - 880 format (/,' Dipole Polarizability Parameters :', + write (iout,910) + 910 format (/,' Dipole Polarizability Parameters :', & //,10x,'Atom Number',5x,'Alpha', & 6x,'Polarization Group',/) end if end if if (use_tholed) then - write (iout,890) i,ia,polarity(ia),thole(ia), + write (iout,920) i,ia,polarity(ia),thole(ia), & tholed(ia),(ip11(j,ia),j=1,np11(ia)) - 890 format (i6,3x,i6,6x,f10.4,2f9.3,3x,120i6) + 920 format (i6,3x,i6,6x,f10.4,2f9.3,3x,120i6) else if (use_thole) then - write (iout,900) i,ia,polarity(ia),thole(ia), + write (iout,930) i,ia,polarity(ia),thole(ia), & (ip11(j,ia),j=1,np11(ia)) - 900 format (i6,3x,i6,6x,f10.4,f9.3,3x,120i6) + 930 format (i6,3x,i6,6x,f10.4,f9.3,3x,120i6) else - write (iout,910) i,ia,polarity(ia), + write (iout,940) i,ia,polarity(ia), & (ip11(j,ia),j=1,np11(ia)) - 910 format (i6,3x,i6,6x,f10.4,3x,120i6) + 940 format (i6,3x,i6,6x,f10.4,3x,120i6) end if end if end do @@ -2085,13 +2120,13 @@ subroutine paramyze (active) if (active(ia)) then if (header) then header = .false. - write (iout,920) - 920 format (/,' Charge Transfer Parameters :', + write (iout,950) + 950 format (/,' Charge Transfer Parameters :', & //,10x,'Atom Number',23x,'Charge', & 6x,'Damp',/) end if - write (iout,930) i,ia,chgct(ia),dmpct(ia) - 930 format (i6,3x,i6,25x,f10.3,f10.4) + write (iout,960) i,ia,chgct(ia),dmpct(ia) + 960 format (i6,3x,i6,25x,f10.3,f10.4) end if end do end if @@ -2108,13 +2143,13 @@ subroutine paramyze (active) if (active(ia) .or. active(ib)) then if (header) then header = .false. - write (iout,940) - 940 format (/,' Bond Charge Flux Parameters :', + write (iout,970) + 970 format (/,' Bond Charge Flux Parameters :', & //,10x,'Atom Numbers',24x,'KCFB',/) end if k = k + 1 - write (iout,950) k,ia,ib,bflx(i) - 950 format (i6,3x,2i6,19x,f10.4) + write (iout,980) k,ia,ib,bflx(i) + 980 format (i6,3x,2i6,19x,f10.4) end if end if end do @@ -2134,15 +2169,15 @@ subroutine paramyze (active) if (active(ia) .or. active(ib) .or. active(ic)) then if (header) then header = .false. - write (iout,960) - 960 format (/,' Angle Charge Flux Parameters :', + write (iout,990) + 990 format (/,' Angle Charge Flux Parameters :', & //,13x,'Atom Numbers',17x,'KACF1', & 5x,'KACF2',5x,'KBCF1',5x,'KBCF2',/) end if k = k + 1 - write (iout,970) k,ia,ib,ic,aflx(1,i),aflx(2,i), + write (iout,1000) k,ia,ib,ic,aflx(1,i),aflx(2,i), & abflx(1,i),abflx(2,i) - 970 format (i6,3x,3i6,10x,4f10.4) + 1000 format (i6,3x,3i6,10x,4f10.4) end if end if end do @@ -2159,34 +2194,34 @@ subroutine paramyze (active) if (header) then header = .false. if (solvtyp(1:2).eq.'GK') then - write (iout,980) - 980 format (/,' Implicit Solvation Parameters :', + write (iout,1010) + 1010 format (/,' Implicit Solvation Parameters :', & //,10x,'Atom Number',11x,'Rsolv', & 4x,'Rdescr',5x,'S-HCT',4x,'S-Neck', & 3x,'Surface',/) else if (solvtyp(1:2).eq.'PB') then - write (iout,990) - 990 format (/,' Implicit Solvation Parameters :', + write (iout,1020) + 1020 format (/,' Implicit Solvation Parameters :', & //,10x,'Atom Number',10x,'Radius', & 5x,'S-HCT',4x,'S-Neck',3x,'Surface',/) else - write (iout,1000) - 1000 format (/,' Implicit Solvation Parameters :', + write (iout,1030) + 1030 format (/,' Implicit Solvation Parameters :', & //,10x,'Atom Number',10x,'Radius', & 3x,'Surface',/) end if end if if (solvtyp(1:2).eq.'GK') then - write (iout,1010) k,i,rsolv(i),rdescr(i),shct(i), + write (iout,1040) k,i,rsolv(i),rdescr(i),shct(i), & sneck(i),asolv(i) - 1010 format (i6,3x,i6,12x,5f10.4) + 1040 format (i6,3x,i6,12x,5f10.4) else if (solvtyp(1:2).eq.'PB') then - write (iout,1020) k,i,rsolv(i),shct(i),sneck(i), + write (iout,1050) k,i,rsolv(i),shct(i),sneck(i), & asolv(i) - 1020 format (i6,3x,i6,12x,4f10.4) + 1050 format (i6,3x,i6,12x,4f10.4) else - write (iout,1030) k,i,rsolv(i),asolv(i) - 1030 format (i6,3x,i6,12x,2f10.4) + write (iout,1060) k,i,rsolv(i),asolv(i) + 1060 format (i6,3x,i6,12x,2f10.4) end if end if end do @@ -2201,13 +2236,13 @@ subroutine paramyze (active) j = class(ia) if (header) then header = .false. - write (iout,1040) - 1040 format (/,' Conjugated Pi-Atom Parameters :', + write (iout,1070) + 1070 format (/,' Conjugated Pi-Atom Parameters :', & //,10x,'Atom Number',14x,'Nelect', & 6x,'Ionize',4x,'Repulsion',/) end if - write (iout,1050) i,ia,electron(j),ionize(j),repulse(j) - 1050 format (i6,3x,i6,17x,f8.1,3x,f10.4,2x,f10.4) + write (iout,1080) i,ia,electron(j),ionize(j),repulse(j) + 1080 format (i6,3x,i6,17x,f8.1,3x,f10.4,2x,f10.4) end do end if c @@ -2220,13 +2255,13 @@ subroutine paramyze (active) ib = ibpi(3,i) if (header) then header = .false. - write (iout,1060) - 1060 format (/,' Conjugated Pi-Bond Parameters :', + write (iout,1090) + 1090 format (/,' Conjugated Pi-Bond Parameters :', & //,10x,'Atom Numbers',21x,'K Slope', & 3x,'L Slope',/) end if - write (iout,1070) i,ia,ib,kslope(i),lslope(i) - 1070 format (i6,3x,2i6,19x,2f10.4) + write (iout,1100) i,ia,ib,kslope(i),lslope(i) + 1100 format (i6,3x,2i6,19x,2f10.4) end do end if return diff --git a/source/chkpole.f b/source/chkpole.f index 31c1e4b2..38f7cf6a 100644 --- a/source/chkpole.f +++ b/source/chkpole.f @@ -20,6 +20,7 @@ subroutine chkpole use atoms use mpole use repel + use xrepel implicit none integer i,k integer ia,ib,ic,id @@ -28,6 +29,7 @@ subroutine chkpole real*8 xcd,ycd,zcd real*8 c1,c2,c3,vol logical dopol,dorep + logical doxrep logical check c c @@ -41,8 +43,10 @@ subroutine chkpole end if if (allocated(replist)) then if (replist(i) .ne. 0) dorep = .true. + else if (allocated(xreplist)) then + if (xreplist(i) .ne. 0) doxrep = .true. end if - if (dopol .or. dorep) then + if (dopol .or. dorep .or. doxrep) then check = .true. if (polaxe(i) .ne. 'Z-then-X') check = .false. if (yaxis(i) .eq. 0) check = .false. @@ -87,6 +91,8 @@ subroutine chkpole repole(8,i) = -repole(8,i) repole(10,i) = -repole(10,i) repole(12,i) = -repole(12,i) + else if (doxrep) then + xrepole(3,i) = -xrepole(3,i) end if end if end if diff --git a/source/energy.f b/source/energy.f index f00c3105..2ce5378e 100644 --- a/source/energy.f +++ b/source/energy.f @@ -116,7 +116,11 @@ function energy () if (vdwtyp .eq. 'BUFFERED-14-7') call ehal if (vdwtyp .eq. 'GAUSSIAN') call egauss end if - if (use_repel) call erepel + if (use_repel) then + call erepel + else if (use_xrepel) then + call exrepel + end if if (use_disp) call edisp c c call any miscellaneous energy component routines diff --git a/source/exrepel.f b/source/exrepel.f new file mode 100644 index 00000000..3995fdd8 --- /dev/null +++ b/source/exrepel.f @@ -0,0 +1,698 @@ +c +c +c ############################################################ +c ## COPYRIGHT (C) 2022 by Moses KJ Chung & Jay W. Ponder ## +c ## All Rights Reserved ## +c ############################################################ +c +c ######################################################### +c ## ## +c ## subroutine exrepel -- exchange repulsion energy ## +c ## ## +c ######################################################### +c +c +c "exrepel" calculates the exchange repulsion energy +c +c literature reference: +c +c M. K. J. Chung and J. W. Ponder, "Beyond isotropic repulsion: +c Classical anisotropic repulsion by inclusion of p orbitals", +c Journal of Chemical Physics, 160, 174118 (2024) +c +c + subroutine exrepel + use limits + implicit none +c +c +c choose the method for summing over pairwise interactions +c + if (use_mlist) then + call exrepel0b + else + call exrepel0a + end if + return + end +c +c +c ################################################################ +c ## ## +c ## subroutine exrepel0a -- exch repulsion energy via loop ## +c ## ## +c ################################################################ +c +c +c "exrepel0a" calculates the exchange repulsion energy using a +c double loop +c +c + subroutine exrepel0a + use atoms + use bound + use cell + use couple + use energi + use group + use mutant + use polar + use reppot + use shunt + use units + use usage + use xrepel + implicit none + integer i,j,k + integer ii,kk + integer ind1,ind2,ind3 + integer jcell + real*8 e,fgrp,taper + real*8 xi,yi,zi + real*8 xr,yr,zr + real*8 r,r2,r3,r4,r5 + real*8 normi + real*8 zxri,zxrk + real*8 vali,valk + real*8 dmpi,dmpk + real*8 dis,dks + real*8 dmpip,dmpkp + real*8 cis,cks + real*8 cix,ckx + real*8 ciy,cky + real*8 ciz,ckz + real*8 rcix,rckx + real*8 rciy,rcky + real*8 rciz,rckz + real*8 cscs,cxcx,cycy + real*8 czcz,cscz,czcs + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 intS,intS2 + real*8 bi(3) + real*8 bj(3) + real*8 bk(3) + real*8, allocatable :: rscale(:) + logical proceed,usei + logical muti,mutk,mutik + logical grad + character*6 mode +c +c +c zero out the repulsion energy contribution +c + er = 0.0d0 + if (nxrep .eq. 0) return +c +c check the sign of multipole components at chiral sites +c + call chkpole +c +c determine pseudo orbital coefficients +c + call solvcoeff +c +c rotate the coefficient components into the global frame +c + call rotcoeff +c +c perform dynamic allocation of some local arrays +c + allocate (rscale(n)) +c +c initialize connected atom exclusion coefficients +c + do i = 1, n + rscale(i) = 1.0d0 + end do +c +c set cutoff distances and switching coefficients +c + mode = 'REPULS' + call switch (mode) +c +c set gradient mode to false +c + grad = .false. +c +c calculate the exchange repulsion interaction energy term +c + do ii = 1, nxrep-1 + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kk = ii+1, nxrep + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + if (use_bounds) call image (xr,yr,zr) + r2 = xr * xr + yr * yr + zr * zr + if (r2 .le. off2) then + r = sqrt(r2) + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, 0.0d0, + & r, grad, SS, dSS, SPz, dSPz, PzS, dPzS, + & PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + e = hartree*(zxri*valk+zxrk*vali)*intS2/r*rscale(k) +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + r3 = r2 * r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + e = e * taper + end if +c +c scale the interaction based on its group membership +c + if (use_group) e = e * fgrp +c +c increment the overall exchange repulsion energy component +c + er = er + e + end if + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do +c +c for periodic boundary conditions with large cutoffs +c neighbors must be found by the replicates method +c + if (use_replica) then +c +c calculate interaction energy with other unit cells +c + do ii = 1, nxrep + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kk = ii, nxrep + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + do jcell = 2, ncell + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + call imager (xr,yr,zr,jcell) + r2 = xr*xr + yr* yr + zr*zr + if (r2 .le. off2) then + r = sqrt(r2) + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, + & 0.0d0, r, grad, SS, dSS, SPz, dSPz, PzS, + & dPzS, PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + e = hartree*(zxri*valk+zxrk*vali)*intS2/r* + & rscale(k) +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + r3 = r2 * r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + e = e * taper + end if +c +c scale the interaction based on its group membership +c + if (use_group) e = e * fgrp +c +c increment the overall exchange repulsion energy component; +c interaction of an atom with its own image counts half +c + if (i .eq. k) e = 0.5d0 * e + er = er + e + end if + end do + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do + end if +c +c perform deallocation of some local arrays +c + deallocate (rscale) + return + end +c +c +c ################################################################ +c ## ## +c ## subroutine exrepel0B -- exch repulsion energy via list ## +c ## ## +c ################################################################ +c +c +c "exrepel0b" calculates the exchange repulsion energy using a +c pairwise neightbor list +c +c + subroutine exrepel0b + use atoms + use bound + use couple + use energi + use group + use mutant + use neigh + use polar + use reppot + use shunt + use units + use usage + use xrepel + implicit none + integer i,j,k + integer ii,kk,kkk + integer ind1,ind2,ind3 + real*8 e,fgrp,taper + real*8 xi,yi,zi + real*8 xr,yr,zr + real*8 r,r2,r3,r4,r5 + real*8 normi + real*8 zxri,zxrk + real*8 vali,valk + real*8 dmpi,dmpk + real*8 dis,dks + real*8 dmpip,dmpkp + real*8 cis,cks + real*8 cix,ckx + real*8 ciy,cky + real*8 ciz,ckz + real*8 rcix,rckx + real*8 rciy,rcky + real*8 rciz,rckz + real*8 cscs,cxcx,cycy + real*8 czcz,cscz,czcs + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 intS,intS2 + real*8 bi(3) + real*8 bj(3) + real*8 bk(3) + real*8, allocatable :: rscale(:) + logical proceed,usei + logical muti,mutk,mutik + logical grad + character*6 mode +c +c +c zero out the repulsion energy contribution +c + er = 0.0d0 + if (nxrep .eq. 0) return +c +c check the sign of multipole components at chiral sites +c + call chkpole +c +c determine pseudo orbital coefficients +c + call solvcoeff +c +c rotate the coefficient components into the global frame +c + call rotcoeff +c +c perform dynamic allocation of some local arrays +c + allocate (rscale(n)) +c +c initialize connected atom exclusion coefficients +c + do i = 1, n + rscale(i) = 1.0d0 + end do +c +c set cutoff distances and switching coefficients +c + mode = 'REPULS' + call switch (mode) +c +c set gradient mode to false +c + grad = .false. +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(nxrep,ixrep,x,y,z,zpxr,dmppxr,rcpxr,n12,i12, +!$OMP& n13,i13,n14,i14,n15,i15,r2scale,r3scale,r4scale,r5scale, +!$OMP& nelst,elst,use,use_group,use_intra,use_bounds,grad, +!$OMP& mut,cut2,off2,c0,c1,c2,c3,c4,c5) +!$OMP& firstprivate(rscale) +!$OMP& shared (er) +!$OMP DO reduction(+:er) schedule(guided) +c +c calculate the exchange repulsion interaction energy term +c + do ii = 1, nxrep + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kkk = 1, nelst(ii) + kk = elst(kkk,ii) + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + if (use_bounds) call image (xr,yr,zr) + r2 = xr * xr + yr * yr + zr * zr + if (r2 .le. off2) then + r = sqrt(r2) + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, 0.0d0, + & r, grad, SS, dSS, SPz, dSPz, PzS, dPzS, + & PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + e = hartree*(zxri*valk+zxrk*vali)*intS2/r*rscale(k) +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + r3 = r2 * r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + e = e * taper + end if +c +c scale the interaction based on its group membership +c + if (use_group) e = e * fgrp +c +c increment the overall exchange repulsion energy component +c + er = er + e + end if + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do +c +c OpenMP directives for the major loop structure +c +!$OMP END DO +!$OMP END PARALLEL +c +c perform deallocation of some local arrays +c + deallocate (rscale) + return + end diff --git a/source/exrepel1.f b/source/exrepel1.f new file mode 100644 index 00000000..8f57860d --- /dev/null +++ b/source/exrepel1.f @@ -0,0 +1,1481 @@ +c +c +c ############################################################ +c ## COPYRIGHT (C) 2022 by Moses KJ Chung & Jay W. Ponder ## +c ## All Rights Reserved ## +c ############################################################ +c +c ############################################################### +c ## ## +c ## subroutine exrepel1 -- exch repulsion energy & derivs ## +c ## ## +c ############################################################### +c +c +c "exrepel1" calculates the exchange repulsion energy and first +c derivatives with respect to Cartesian coordinates +c +c literature reference: +c +c M. K. J. Chung and J. W. Ponder, "Beyond isotropic repulsion: +c Classical anisotropic repulsion by inclusion of p orbitals", +c Journal of Chemical Physics, 160, 174118 (2024) +c +c + subroutine exrepel1 + use limits + implicit none +c +c +c choose the method for summing over pairwise interactions +c + if (use_mlist) then + call exrepel1b + else + call exrepel1a + end if + return + end +c +c +c ################################################################## +c ## ## +c ## subroutine exrepel1a -- exch repulsion analysis via loop ## +c ## ## +c ################################################################## +c +c +c "exrepel1a" calculates the exchange repulsion energy and first +c derivatives with respect to Cartesian coordinates using a +c pairwise double loop +c +c + subroutine exrepel1a + use atoms + use bound + use cell + use couple + use deriv + use energi + use group + use mpole + use mutant + use reppot + use shunt + use units + use usage + use virial + use xrepel + implicit none + integer i,j,k + integer ii,kk,jcell + integer ix,iy,iz + integer ind1,ind2,ind3 + real*8 e,fgrp + real*8 taper,dtaper + real*8 xi,yi,zi + real*8 xr,yr,zr + real*8 xix,yix,zix + real*8 xiy,yiy,ziy + real*8 xiz,yiz,ziz + real*8 rr1 + real*8 r,r2,r3,r4,r5 + real*8 normi + real*8 zxri,zxrk + real*8 vali,valk + real*8 dmpi,dmpk + real*8 dis,dks + real*8 dmpip,dmpkp + real*8 cis,cks + real*8 cix,ckx + real*8 ciy,cky + real*8 ciz,ckz + real*8 rcix,rckx + real*8 rciy,rcky + real*8 rciz,rckz + real*8 rcixr,rckxr + real*8 rciyr,rckyr + real*8 rcizr,rckzr + real*8 cscs,cxcx,cycy + real*8 czcz,cscz,czcs + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 intS,intS2 + real*8 dintS + real*8 dintSx,dintSy,dintSz + real*8 intSR,preintSR + real*8 pre + real*8 term1 + real*8 term2x,term2y,term2z + real*8 frcx,frcy,frcz + real*8 ncix,nciy,nciz + real*8 nckx,ncky,nckz + real*8 nrcix,nrciy,nrciz + real*8 nrckx,nrcky,nrckz + real*8 drcixdx,drcixdy,drcixdz + real*8 drciydx,drciydy,drciydz + real*8 drcizdx,drcizdy,drcizdz + real*8 drckxdx,drckxdy,drckxdz + real*8 drckydx,drckydy,drckydz + real*8 drckzdx,drckzdy,drckzdz + real*8 tixintS,tiyintS,tizintS + real*8 tkxintS,tkyintS,tkzintS + real*8 vxx,vyy,vzz + real*8 vxy,vxz,vyz + real*8 bi(3) + real*8 bj(3) + real*8 bk(3) + real*8 ttri(3),ttrk(3) + real*8 fix(3),fiy(3),fiz(3) + real*8, allocatable :: rscale(:) + real*8, allocatable :: ter(:,:) + logical proceed,usei + logical muti,mutk,mutik + logical header,huge + logical grad + character*6 mode +c +c +c zero out the repulsion energy and derivatives +c + er = 0.0d0 + do i = 1, n + der(1,i) = 0.0d0 + der(2,i) = 0.0d0 + der(3,i) = 0.0d0 + end do + if (nxrep .eq. 0) return +c +c check the sign of multipole components at chiral sites +c + call chkpole +c +c determine pseudo orbital coefficients +c + call solvcoeff +c +c rotate the coefficient components into the global frame +c + call rotcoeff +c +c perform dynamic allocation of some local arrays +c + allocate (rscale(n)) + allocate (ter(3,n)) +c +c initialize connected atom exclusion coefficients +c + do i = 1, n + rscale(i) = 1.0d0 + do j = 1, 3 + ter(j,i) = 0.0d0 + end do + end do +c +c set cutoff distances and switching coefficients +c + mode = 'REPULS' + call switch (mode) +c +c set gradient mode to true +c + grad = .true. +c +c calculate the exchange repulsion energy and derivative +c + do ii = 1, nxrep-1 + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kk = ii+1, nxrep + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + if (use_bounds) call image (xr,yr,zr) + r2 = xr * xr + yr * yr + zr * zr + if (r2 .le. off2) then + r = sqrt(r2) + r3 = r2 * r + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, 0.0d0, + & r, grad, SS, dSS, SPz, dSPz, PzS, dPzS, + & PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + dintS = cscs * dSS + cxcx * dPxPx + cycy * dPyPy + & + czcz * dPzPz + cscz * dSPz + czcs * dPzS + dintSx = dintS * bk(1) + dintSy = dintS * bk(2) + dintSz = dintS * bk(3) + rcixr = rcix/r + rciyr = rciy/r + rcizr = rciz/r + rckxr = rckx/r + rckyr = rcky/r + rckzr = rckz/r + drcixdx = bi(1)*(-rcizr) + drcixdy = bi(2)*(-rcizr) + drcixdz = bi(3)*(-rcizr) + drciydx = bj(1)*(-rcizr) + drciydy = bj(2)*(-rcizr) + drciydz = bj(3)*(-rcizr) + drcizdx = bi(1)*( rcixr) + bj(1)*( rciyr) + drcizdy = bi(2)*( rcixr) + bj(2)*( rciyr) + drcizdz = bi(3)*( rcixr) + bj(3)*( rciyr) + drckxdx = bi(1)*(-rckzr) + drckxdy = bi(2)*(-rckzr) + drckxdz = bi(3)*(-rckzr) + drckydx = bj(1)*(-rckzr) + drckydy = bj(2)*(-rckzr) + drckydz = bj(3)*(-rckzr) + drckzdx = bi(1)*( rckxr) + bj(1)*( rckyr) + drckzdy = bi(2)*( rckxr) + bj(2)*( rckyr) + drckzdz = bi(3)*( rckxr) + bj(3)*( rckyr) + dintSx = dintSx + drcizdx*cks*pzs + drcixdx*rckx*pxpx + & + drciydx*rcky*pypy + drcizdx*rckz*pzpz + dintSy = dintSy + drcizdy*cks*pzs + drcixdy*rckx*pxpx + & + drciydy*rcky*pypy + drcizdy*rckz*pzpz + dintSz = dintSz + drcizdz*cks*pzs + drcixdz*rckx*pxpx + & + drciydz*rcky*pypy + drcizdz*rckz*pzpz + dintSx = dintSx + cis*drckzdx*spz + rcix*drckxdx*pxpx + & + rciy*drckydx*pypy + rciz*drckzdx*pzpz + dintSy = dintSy + cis*drckzdy*spz + rcix*drckxdy*pxpx + & + rciy*drckydy*pypy + rciz*drckzdy*pzpz + dintSz = dintSz + cis*drckzdz*spz + rcix*drckxdz*pxpx + & + rciy*drckydz*pypy + rciz*drckzdz*pzpz + pre = hartree * (zxri*valk + zxrk*vali) * rscale(k) + e = pre * intS2 / r + term1 = -intS2 / r3 + intSR = 2.0d0 * intS / r + term2x = intSR * dintSx + term2y = intSR * dintSy + term2z = intSR * dintSz + +c +c compute the force components for this interaction +c + frcx = -pre * (xr * term1 + term2x) + frcy = -pre * (yr * term1 + term2y) + frcz = -pre * (zr * term1 + term2z) +c +c compute the torque components for this interaction +c + ncix = 0.0d0 + nciy = -ciz + nciz = ciy + nrcix = bi(1)*ncix + bi(2)*nciy + bi(3)*nciz + nrciy = bj(1)*ncix + bj(2)*nciy + bj(3)*nciz + nrciz = bk(1)*ncix + bk(2)*nciy + bk(3)*nciz + cscs = 0.0d0 * cks + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + cscz = 0.0d0 * rckz + czcs = nrciz * cks + tixintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + ncix = ciz + nciy = 0.0d0 + nciz = -cix + nrcix = bi(1)*ncix + bi(2)*nciy + bi(3)*nciz + nrciy = bj(1)*ncix + bj(2)*nciy + bj(3)*nciz + nrciz = bk(1)*ncix + bk(2)*nciy + bk(3)*nciz + cscs = 0.0d0 * cks + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + cscz = 0.0d0 * rckz + czcs = nrciz * cks + tiyintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + ncix = -ciy + nciy = cix + nciz = 0.0d0 + nrcix = bi(1)*ncix + bi(2)*nciy + bi(3)*nciz + nrciy = bj(1)*ncix + bj(2)*nciy + bj(3)*nciz + nrciz = bk(1)*ncix + bk(2)*nciy + bk(3)*nciz + cscs = 0.0d0 * cks + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + cscz = 0.0d0 * rckz + czcs = nrciz * cks + tizintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + nckx = 0.0d0 + ncky = -ckz + nckz = cky + nrckx = bi(1)*nckx + bi(2)*ncky + bi(3)*nckz + nrcky = bj(1)*nckx + bj(2)*ncky + bj(3)*nckz + nrckz = bk(1)*nckx + bk(2)*ncky + bk(3)*nckz + cscs = cis * 0.0d0 + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + czcs = rciz * 0.0d0 + tkxintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + nckx = ckz + ncky = 0.0d0 + nckz = -ckx + nrckx = bi(1)*nckx + bi(2)*ncky + bi(3)*nckz + nrcky = bj(1)*nckx + bj(2)*ncky + bj(3)*nckz + nrckz = bk(1)*nckx + bk(2)*ncky + bk(3)*nckz + cscs = cis * 0.0d0 + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + czcs = rciz * 0.0d0 + tkyintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + nckx = -cky + ncky = ckx + nckz = 0.0d0 + nrckx = bi(1)*nckx + bi(2)*ncky + bi(3)*nckz + nrcky = bj(1)*nckx + bj(2)*ncky + bj(3)*nckz + nrckz = bk(1)*nckx + bk(2)*ncky + bk(3)*nckz + cscs = cis * 0.0d0 + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + czcs = rciz * 0.0d0 + tkzintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + preintSR = -pre * intSR + ttri(1) = preintSR * tixintS + ttri(2) = preintSR * tiyintS + ttri(3) = preintSR * tizintS + ttrk(1) = preintSR * tkxintS + ttrk(2) = preintSR * tkyintS + ttrk(3) = preintSR * tkzintS +c +c scale the interaction based on its group membership +c + if (use_group) then + e = fgrp * e + frcx = fgrp * frcx + frcy = fgrp * frcy + frcz = fgrp * frcz + do j = 1, 3 + ttri(j) = fgrp * ttri(j) + ttrk(j) = fgrp * ttrk(j) + end do + end if +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + rr1 = 1.0d0 / r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 + & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 + dtaper = dtaper * e * rr1 + e = e * taper + frcx = frcx*taper - dtaper*xr + frcy = frcy*taper - dtaper*yr + frcz = frcz*taper - dtaper*zr + do j = 1, 3 + ttri(j) = ttri(j) * taper + ttrk(j) = ttrk(j) * taper + end do + end if +c +c increment the overall exchange repulsion energy component +c + er = er + e +c +c increment force-based gradient on first site +c + der(1,i) = der(1,i) + frcx + der(2,i) = der(2,i) + frcy + der(3,i) = der(3,i) + frcz + ter(1,i) = ter(1,i) + ttri(1) + ter(2,i) = ter(2,i) + ttri(2) + ter(3,i) = ter(3,i) + ttri(3) +c +c increment force-based gradient and torque on second site +c + der(1,k) = der(1,k) - frcx + der(2,k) = der(2,k) - frcy + der(3,k) = der(3,k) - frcz + ter(1,k) = ter(1,k) + ttrk(1) + ter(2,k) = ter(2,k) + ttrk(2) + ter(3,k) = ter(3,k) + ttrk(3) +c +c increment the virial due to pairwise Cartesian forces +c + vxx = -xr * frcx + vxy = -0.5d0 * (yr*frcx+xr*frcy) + vxz = -0.5d0 * (zr*frcx+xr*frcz) + vyy = -yr * frcy + vyz = -0.5d0 * (zr*frcy+yr*frcz) + vzz = -zr * frcz + vir(1,1) = vir(1,1) + vxx + vir(2,1) = vir(2,1) + vxy + vir(3,1) = vir(3,1) + vxz + vir(1,2) = vir(1,2) + vxy + vir(2,2) = vir(2,2) + vyy + vir(3,2) = vir(3,2) + vyz + vir(1,3) = vir(1,3) + vxz + vir(2,3) = vir(2,3) + vyz + vir(3,3) = vir(3,3) + vzz + end if + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do +c +c for periodic boundary conditions with large cutoffs +c neighbors must be found by the replicates method +c + if (use_replica) then +c +c calculate interaction energy with other unit cells +c + do ii = 1, nxrep + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kk = ii, nxrep + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + do jcell = 2, ncell + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + call imager (xr,yr,zr,jcell) + r2 = xr*xr + yr* yr + zr*zr + if (r2 .le. off2) then + r = sqrt(r2) + r3 = r2 * r + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, + & 0.0d0, r, grad, SS, dSS, SPz, dSPz, PzS, + & dPzS, PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + dintS = cscs * dSS + cxcx * dPxPx + cycy * dPyPy + & + czcz * dPzPz + cscz * dSPz + czcs * dPzS + dintSx = dintS * bk(1) + dintSy = dintS * bk(2) + dintSz = dintS * bk(3) + rcixr = rcix/r + rciyr = rciy/r + rcizr = rciz/r + rckxr = rckx/r + rckyr = rcky/r + rckzr = rckz/r + drcixdx = bi(1)*(-rcizr) + drcixdy = bi(2)*(-rcizr) + drcixdz = bi(3)*(-rcizr) + drciydx = bj(1)*(-rcizr) + drciydy = bj(2)*(-rcizr) + drciydz = bj(3)*(-rcizr) + drcizdx = bi(1)*( rcixr) + bj(1)*( rciyr) + drcizdy = bi(2)*( rcixr) + bj(2)*( rciyr) + drcizdz = bi(3)*( rcixr) + bj(3)*( rciyr) + drckxdx = bi(1)*(-rckzr) + drckxdy = bi(2)*(-rckzr) + drckxdz = bi(3)*(-rckzr) + drckydx = bj(1)*(-rckzr) + drckydy = bj(2)*(-rckzr) + drckydz = bj(3)*(-rckzr) + drckzdx = bi(1)*( rckxr) + bj(1)*( rckyr) + drckzdy = bi(2)*( rckxr) + bj(2)*( rckyr) + drckzdz = bi(3)*( rckxr) + bj(3)*( rckyr) + dintSx =dintSx+drcizdx*cks*pzs+drcixdx*rckx*pxpx + & + drciydx*rcky*pypy + drcizdx*rckz*pzpz + dintSy =dintSy+drcizdy*cks*pzs+drcixdy*rckx*pxpx + & + drciydy*rcky*pypy + drcizdy*rckz*pzpz + dintSz =dintSz+drcizdz*cks*pzs+drcixdz*rckx*pxpx + & + drciydz*rcky*pypy + drcizdz*rckz*pzpz + dintSx =dintSx+cis*drckzdx*spz+rcix*drckxdx*pxpx + & + rciy*drckydx*pypy + rciz*drckzdx*pzpz + dintSy =dintSy+cis*drckzdy*spz+rcix*drckxdy*pxpx + & + rciy*drckydy*pypy + rciz*drckzdy*pzpz + dintSz =dintSz+cis*drckzdz*spz+rcix*drckxdz*pxpx + & + rciy*drckydz*pypy + rciz*drckzdz*pzpz + pre = hartree * (zxri*valk+zxrk*vali)*rscale(k) + e = pre * intS2 / r + term1 = -intS2 / r3 + intSR = 2.0d0 * intS / r + term2x = intSR * dintSx + term2y = intSR * dintSy + term2z = intSR * dintSz +c +c compute the force components for this interaction +c + frcx = -pre * (xr * term1 + term2x) + frcy = -pre * (yr * term1 + term2y) + frcz = -pre * (zr * term1 + term2z) +c +c compute the torque components for this interaction +c + ncix = 0.0d0 + nciy = -ciz + nciz = ciy + nrcix = bi(1)*ncix + bi(2)*nciy + bi(3)*nciz + nrciy = bj(1)*ncix + bj(2)*nciy + bj(3)*nciz + nrciz = bk(1)*ncix + bk(2)*nciy + bk(3)*nciz + cscs = 0.0d0 * cks + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + cscz = 0.0d0 * rckz + czcs = nrciz * cks + tixintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + ncix = ciz + nciy = 0.0d0 + nciz = -cix + nrcix = bi(1)*ncix + bi(2)*nciy + bi(3)*nciz + nrciy = bj(1)*ncix + bj(2)*nciy + bj(3)*nciz + nrciz = bk(1)*ncix + bk(2)*nciy + bk(3)*nciz + cscs = 0.0d0 * cks + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + cscz = 0.0d0 * rckz + czcs = nrciz * cks + tiyintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + ncix = -ciy + nciy = cix + nciz = 0.0d0 + nrcix = bi(1)*ncix + bi(2)*nciy + bi(3)*nciz + nrciy = bj(1)*ncix + bj(2)*nciy + bj(3)*nciz + nrciz = bk(1)*ncix + bk(2)*nciy + bk(3)*nciz + cscs = 0.0d0 * cks + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + cscz = 0.0d0 * rckz + czcs = nrciz * cks + tizintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + nckx = 0.0d0 + ncky = -ckz + nckz = cky + nrckx = bi(1)*nckx + bi(2)*ncky + bi(3)*nckz + nrcky = bj(1)*nckx + bj(2)*ncky + bj(3)*nckz + nrckz = bk(1)*nckx + bk(2)*ncky + bk(3)*nckz + cscs = cis * 0.0d0 + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + czcs = rciz * 0.0d0 + tkxintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + nckx = ckz + ncky = 0.0d0 + nckz = -ckx + nrckx = bi(1)*nckx + bi(2)*ncky + bi(3)*nckz + nrcky = bj(1)*nckx + bj(2)*ncky + bj(3)*nckz + nrckz = bk(1)*nckx + bk(2)*ncky + bk(3)*nckz + cscs = cis * 0.0d0 + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + czcs = rciz * 0.0d0 + tkyintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + nckx = -cky + ncky = ckx + nckz = 0.0d0 + nrckx = bi(1)*nckx + bi(2)*ncky + bi(3)*nckz + nrcky = bj(1)*nckx + bj(2)*ncky + bj(3)*nckz + nrckz = bk(1)*nckx + bk(2)*ncky + bk(3)*nckz + cscs = cis * 0.0d0 + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + czcs = rciz * 0.0d0 + tkzintS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + preintSR = -pre * intSR + ttri(1) = preintSR * tixintS + ttri(2) = preintSR * tiyintS + ttri(3) = preintSR * tizintS + ttrk(1) = preintSR * tkxintS + ttrk(2) = preintSR * tkyintS + ttrk(3) = preintSR * tkzintS +c +c scale the interaction based on its group membership +c + if (use_group) then + e = fgrp * e + frcx = fgrp * frcx + frcy = fgrp * frcy + frcz = fgrp * frcz + do j = 1, 3 + ttri(j) = fgrp * ttri(j) + ttrk(j) = fgrp * ttrk(j) + end do + end if +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + rr1 = 1.0d0 / r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 + & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 + dtaper = dtaper * e * rr1 + e = e * taper + frcx = frcx*taper - dtaper*xr + frcy = frcy*taper - dtaper*yr + frcz = frcz*taper - dtaper*zr + do j = 1, 3 + ttri(j) = ttri(j) * taper + ttrk(j) = ttrk(j) * taper + end do + end if +c +c increment the overall exchange repulsion energy component; +c interaction of an atom with its own image counts half +c + if (i .eq. k) then + e = 0.5d0 * e + frcx = 0.5d0 * frcx + frcy = 0.5d0 * frcy + frcz = 0.5d0 * frcz + do j = 1, 3 + ttri(j) = 0.5d0 * ttri(j) + ttrk(j) = 0.5d0 * ttrk(j) + end do + end if +c +c increment the overall exchange repulsion energy component +c + er = er + e +c +c increment force-based gradient on first site +c + der(1,i) = der(1,i) + frcx + der(2,i) = der(2,i) + frcy + der(3,i) = der(3,i) + frcz + ter(1,i) = ter(1,i) + ttri(1) + ter(2,i) = ter(2,i) + ttri(2) + ter(3,i) = ter(3,i) + ttri(3) +c +c increment force-based gradient and torque on second site +c + der(1,k) = der(1,k) - frcx + der(2,k) = der(2,k) - frcy + der(3,k) = der(3,k) - frcz + ter(1,k) = ter(1,k) + ttrk(1) + ter(2,k) = ter(2,k) + ttrk(2) + ter(3,k) = ter(3,k) + ttrk(3) +c +c increment the virial due to pairwise Cartesian forces +c + vxx = -xr * frcx + vxy = -0.5d0 * (yr*frcx+xr*frcy) + vxz = -0.5d0 * (zr*frcx+xr*frcz) + vyy = -yr * frcy + vyz = -0.5d0 * (zr*frcy+yr*frcz) + vzz = -zr * frcz + vir(1,1) = vir(1,1) + vxx + vir(2,1) = vir(2,1) + vxy + vir(3,1) = vir(3,1) + vxz + vir(1,2) = vir(1,2) + vxy + vir(2,2) = vir(2,2) + vyy + vir(3,2) = vir(3,2) + vyz + vir(1,3) = vir(1,3) + vxz + vir(2,3) = vir(2,3) + vyz + vir(3,3) = vir(3,3) + vzz + end if + end do + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do + end if +c +c resolve site torques then increment forces and virial +c + do ii = 1, nxrep + i = ixrep(ii) + call torque (i,ter(1,i),fix,fiy,fiz,der) + iz = zaxis(i) + ix = xaxis(i) + iy = abs(yaxis(i)) + if (iz .eq. 0) iz = i + if (ix .eq. 0) ix = i + if (iy .eq. 0) iy = i + xiz = x(iz) - x(i) + yiz = y(iz) - y(i) + ziz = z(iz) - z(i) + xix = x(ix) - x(i) + yix = y(ix) - y(i) + zix = z(ix) - z(i) + xiy = x(iy) - x(i) + yiy = y(iy) - y(i) + ziy = z(iy) - z(i) + vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) + vxy = 0.5d0 * (yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) + & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) + vxz = 0.5d0 * (zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) + & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) + vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) + vyz = 0.5d0 * (zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) + & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) + vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) + vir(1,1) = vir(1,1) + vxx + vir(2,1) = vir(2,1) + vxy + vir(3,1) = vir(3,1) + vxz + vir(1,2) = vir(1,2) + vxy + vir(2,2) = vir(2,2) + vyy + vir(3,2) = vir(3,2) + vyz + vir(1,3) = vir(1,3) + vxz + vir(2,3) = vir(2,3) + vyz + vir(3,3) = vir(3,3) + vzz + end do +c +c perform deallocation of some local arrays +c + deallocate (rscale) + deallocate (ter) + return + end +c +c +c ################################################################## +c ## ## +c ## subroutine exrepel1b -- exch repulsion analysis via list ## +c ## ## +c ################################################################## +c +c +c "exrepel1b" calculates the exchange repulsion energy and first +c derivatives with respect to Cartesian coordinates using a +c pairwise neightbor list +c +c + subroutine exrepel1b + use atoms + use bound + use couple + use deriv + use energi + use group + use mpole + use mutant + use neigh + use reppot + use shunt + use units + use usage + use virial + use xrepel + implicit none + integer i,j,k + integer ii,kk,kkk + integer ix,iy,iz + integer ind1,ind2,ind3 + real*8 e,fgrp + real*8 taper,dtaper + real*8 xi,yi,zi + real*8 xr,yr,zr + real*8 xix,yix,zix + real*8 xiy,yiy,ziy + real*8 xiz,yiz,ziz + real*8 rr1 + real*8 r,r2,r3,r4,r5 + real*8 normi + real*8 zxri,zxrk + real*8 vali,valk + real*8 dmpi,dmpk + real*8 dis,dks + real*8 dmpip,dmpkp + real*8 cis,cks + real*8 cix,ckx + real*8 ciy,cky + real*8 ciz,ckz + real*8 rcix,rckx + real*8 rciy,rcky + real*8 rciz,rckz + real*8 rcixr,rckxr + real*8 rciyr,rckyr + real*8 rcizr,rckzr + real*8 cscs,cxcx,cycy + real*8 czcz,cscz,czcs + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 intS,intS2 + real*8 dintS + real*8 dintSx,dintSy,dintSz + real*8 intSR,preintSR + real*8 pre + real*8 term1 + real*8 term2x,term2y,term2z + real*8 frcx,frcy,frcz + real*8 ncix,nciy,nciz + real*8 nckx,ncky,nckz + real*8 nrcix,nrciy,nrciz + real*8 nrckx,nrcky,nrckz + real*8 drcixdx,drcixdy,drcixdz + real*8 drciydx,drciydy,drciydz + real*8 drcizdx,drcizdy,drcizdz + real*8 drckxdx,drckxdy,drckxdz + real*8 drckydx,drckydy,drckydz + real*8 drckzdx,drckzdy,drckzdz + real*8 tixintS,tiyintS,tizintS + real*8 tkxintS,tkyintS,tkzintS + real*8 vxx,vyy,vzz + real*8 vxy,vxz,vyz + real*8 bi(3) + real*8 bj(3) + real*8 bk(3) + real*8 ttri(3),ttrk(3) + real*8 fix(3),fiy(3),fiz(3) + real*8, allocatable :: rscale(:) + real*8, allocatable :: ter(:,:) + logical proceed,usei + logical muti,mutk,mutik + logical header,huge + logical grad + character*6 mode +c +c +c zero out the repulsion energy and derivatives +c + er = 0.0d0 + do i = 1, n + der(1,i) = 0.0d0 + der(2,i) = 0.0d0 + der(3,i) = 0.0d0 + end do + if (nxrep .eq. 0) return +c +c check the sign of multipole components at chiral sites +c + call chkpole +c +c determine pseudo orbital coefficients +c + call solvcoeff +c +c rotate the coefficient components into the global frame +c + call rotcoeff +c +c perform dynamic allocation of some local arrays +c + allocate (rscale(n)) + allocate (ter(3,n)) +c +c initialize connected atom exclusion coefficients +c + do i = 1, n + rscale(i) = 1.0d0 + do j = 1, 3 + ter(j,i) = 0.0d0 + end do + end do +c +c set cutoff distances and switching coefficients +c + mode = 'REPULS' + call switch (mode) +c +c set gradient mode to true +c + grad = .true. +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(nxrep,ixrep,x,y,z,zpxr,dmppxr,rcpxr,n12,i12, +!$OMP& n13,i13,n14,i14,n15,i15,r2scale,r3scale,r4scale,r5scale, +!$OMP& nelst,elst,use,use_group,use_intra,use_bounds,grad, +!$OMP& mut,cut2,off2,xaxis,yaxis,zaxis, +!$OMP& c0,c1,c2,c3,c4,c5) +!$OMP& firstprivate(rscale) shared (er,der,ter,vir) +!$OMP DO reduction(+:er,der,ter,vir) schedule(guided) +c +c calculate the exchange repulsion energy and derivatives +c + do ii = 1, nxrep + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kkk = 1, nelst(ii) + kk = elst(kkk,ii) + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + if (use_bounds) call image (xr,yr,zr) + r2 = xr * xr + yr * yr + zr * zr + if (r2 .le. off2) then + r = sqrt(r2) + r3 = r2 * r + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, 0.0d0, + & r, grad, SS, dSS, SPz, dSPz, PzS, dPzS, + & PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + dintS = cscs * dSS + cxcx * dPxPx + cycy * dPyPy + & + czcz * dPzPz + cscz * dSPz + czcs * dPzS + dintSx = dintS * bk(1) + dintSy = dintS * bk(2) + dintSz = dintS * bk(3) + rcixr = rcix/r + rciyr = rciy/r + rcizr = rciz/r + rckxr = rckx/r + rckyr = rcky/r + rckzr = rckz/r + drcixdx = bi(1)*(-rcizr) + drcixdy = bi(2)*(-rcizr) + drcixdz = bi(3)*(-rcizr) + drciydx = bj(1)*(-rcizr) + drciydy = bj(2)*(-rcizr) + drciydz = bj(3)*(-rcizr) + drcizdx = bi(1)*( rcixr) + bj(1)*( rciyr) + drcizdy = bi(2)*( rcixr) + bj(2)*( rciyr) + drcizdz = bi(3)*( rcixr) + bj(3)*( rciyr) + drckxdx = bi(1)*(-rckzr) + drckxdy = bi(2)*(-rckzr) + drckxdz = bi(3)*(-rckzr) + drckydx = bj(1)*(-rckzr) + drckydy = bj(2)*(-rckzr) + drckydz = bj(3)*(-rckzr) + drckzdx = bi(1)*( rckxr) + bj(1)*( rckyr) + drckzdy = bi(2)*( rckxr) + bj(2)*( rckyr) + drckzdz = bi(3)*( rckxr) + bj(3)*( rckyr) + dintSx = dintSx + drcizdx*cks*pzs + drcixdx*rckx*pxpx + & + drciydx*rcky*pypy + drcizdx*rckz*pzpz + dintSy = dintSy + drcizdy*cks*pzs + drcixdy*rckx*pxpx + & + drciydy*rcky*pypy + drcizdy*rckz*pzpz + dintSz = dintSz + drcizdz*cks*pzs + drcixdz*rckx*pxpx + & + drciydz*rcky*pypy + drcizdz*rckz*pzpz + dintSx = dintSx + cis*drckzdx*spz + rcix*drckxdx*pxpx + & + rciy*drckydx*pypy + rciz*drckzdx*pzpz + dintSy = dintSy + cis*drckzdy*spz + rcix*drckxdy*pxpx + & + rciy*drckydy*pypy + rciz*drckzdy*pzpz + dintSz = dintSz + cis*drckzdz*spz + rcix*drckxdz*pxpx + & + rciy*drckydz*pypy + rciz*drckzdz*pzpz + pre = hartree * (zxri*valk + zxrk*vali) * rscale(k) + e = pre * intS2 / r + term1 = -intS2 / r3 + intSR = 2.0d0 * intS / r + term2x = intSR * dintSx + term2y = intSR * dintSy + term2z = intSR * dintSz + +c +c compute the force components for this interaction +c + frcx = -pre * (xr * term1 + term2x) + frcy = -pre * (yr * term1 + term2y) + frcz = -pre * (zr * term1 + term2z) +c +c compute the torque components for this interaction +c + nciy = -ciz + nciz = ciy + nrcix = bi(2)*nciy + bi(3)*nciz + nrciy = bj(2)*nciy + bj(3)*nciz + nrciz = bk(2)*nciy + bk(3)*nciz + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + czcs = nrciz * cks + tixintS = cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + czcs * PzS + ncix = ciz + nciz = -cix + nrcix = bi(1)*ncix + bi(3)*nciz + nrciy = bj(1)*ncix + bj(3)*nciz + nrciz = bk(1)*ncix + bk(3)*nciz + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + czcs = nrciz * cks + tiyintS = cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + czcs * PzS + ncix = -ciy + nciy = cix + nrcix = bi(1)*ncix + bi(2)*nciy + nrciy = bj(1)*ncix + bj(2)*nciy + nrciz = bk(1)*ncix + bk(2)*nciy + cxcx = nrcix * rckx + cycy = nrciy * rcky + czcz = nrciz * rckz + czcs = nrciz * cks + tizintS = cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + czcs * PzS + ncky = -ckz + nckz = cky + nrckx = bi(2)*ncky + bi(3)*nckz + nrcky = bj(2)*ncky + bj(3)*nckz + nrckz = bk(2)*ncky + bk(3)*nckz + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + tkxintS = cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + nckx = ckz + nckz = -ckx + nrckx = bi(1)*nckx + bi(3)*nckz + nrcky = bj(1)*nckx + bj(3)*nckz + nrckz = bk(1)*nckx + bk(3)*nckz + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + tkyintS = cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + nckx = -cky + ncky = ckx + nrckx = bi(1)*nckx + bi(2)*ncky + nrcky = bj(1)*nckx + bj(2)*ncky + nrckz = bk(1)*nckx + bk(2)*ncky + cxcx = rcix * nrckx + cycy = rciy * nrcky + czcz = rciz * nrckz + cscz = cis * nrckz + tkzintS = cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + preintSR = -pre * intSR + ttri(1) = preintSR * tixintS + ttri(2) = preintSR * tiyintS + ttri(3) = preintSR * tizintS + ttrk(1) = preintSR * tkxintS + ttrk(2) = preintSR * tkyintS + ttrk(3) = preintSR * tkzintS +c +c scale the interaction based on its group membership +c + if (use_group) then + e = fgrp * e + frcx = fgrp * frcx + frcy = fgrp * frcy + frcz = fgrp * frcz + do j = 1, 3 + ttri(j) = fgrp * ttri(j) + ttrk(j) = fgrp * ttrk(j) + end do + end if +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + rr1 = 1.0d0 / r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 + & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 + dtaper = dtaper * e * rr1 + e = e * taper + frcx = frcx*taper - dtaper*xr + frcy = frcy*taper - dtaper*yr + frcz = frcz*taper - dtaper*zr + do j = 1, 3 + ttri(j) = ttri(j) * taper + ttrk(j) = ttrk(j) * taper + end do + end if +c +c increment the overall exchange repulsion energy component +c + er = er + e +c +c increment force-based gradient on first site +c + der(1,i) = der(1,i) + frcx + der(2,i) = der(2,i) + frcy + der(3,i) = der(3,i) + frcz + ter(1,i) = ter(1,i) + ttri(1) + ter(2,i) = ter(2,i) + ttri(2) + ter(3,i) = ter(3,i) + ttri(3) +c +c increment force-based gradient and torque on second site +c + der(1,k) = der(1,k) - frcx + der(2,k) = der(2,k) - frcy + der(3,k) = der(3,k) - frcz + ter(1,k) = ter(1,k) + ttrk(1) + ter(2,k) = ter(2,k) + ttrk(2) + ter(3,k) = ter(3,k) + ttrk(3) +c +c increment the virial due to pairwise Cartesian forces +c + vxx = -xr * frcx + vxy = -0.5d0 * (yr*frcx+xr*frcy) + vxz = -0.5d0 * (zr*frcx+xr*frcz) + vyy = -yr * frcy + vyz = -0.5d0 * (zr*frcy+yr*frcz) + vzz = -zr * frcz + vir(1,1) = vir(1,1) + vxx + vir(2,1) = vir(2,1) + vxy + vir(3,1) = vir(3,1) + vxz + vir(1,2) = vir(1,2) + vxy + vir(2,2) = vir(2,2) + vyy + vir(3,2) = vir(3,2) + vyz + vir(1,3) = vir(1,3) + vxz + vir(2,3) = vir(2,3) + vyz + vir(3,3) = vir(3,3) + vzz + end if + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do +c +c OpenMP directives for the major loop structure +c +!$OMP END DO +!$OMP DO reduction(+:der,vir) schedule(guided) +c +c resolve site torques then increment forces and virial +c + do ii = 1, nxrep + i = ixrep(ii) + call torque (i,ter(1,i),fix,fiy,fiz,der) + iz = zaxis(i) + ix = xaxis(i) + iy = abs(yaxis(i)) + if (iz .eq. 0) iz = i + if (ix .eq. 0) ix = i + if (iy .eq. 0) iy = i + xiz = x(iz) - x(i) + yiz = y(iz) - y(i) + ziz = z(iz) - z(i) + xix = x(ix) - x(i) + yix = y(ix) - y(i) + zix = z(ix) - z(i) + xiy = x(iy) - x(i) + yiy = y(iy) - y(i) + ziy = z(iy) - z(i) + vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) + vxy = 0.5d0 * (yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) + & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) + vxz = 0.5d0 * (zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) + & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) + vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) + vyz = 0.5d0 * (zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) + & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) + vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) + vir(1,1) = vir(1,1) + vxx + vir(2,1) = vir(2,1) + vxy + vir(3,1) = vir(3,1) + vxz + vir(1,2) = vir(1,2) + vxy + vir(2,2) = vir(2,2) + vyy + vir(3,2) = vir(3,2) + vyz + vir(1,3) = vir(1,3) + vxz + vir(2,3) = vir(2,3) + vyz + vir(3,3) = vir(3,3) + vzz + end do +c +c OpenMP directives for the major loop structure +c +!$OMP END DO +!$OMP END PARALLEL +c +c perform deallocation of some local arrays +c + deallocate (rscale) + deallocate (ter) + return + end diff --git a/source/exrepel3.f b/source/exrepel3.f new file mode 100644 index 00000000..03927960 --- /dev/null +++ b/source/exrepel3.f @@ -0,0 +1,1469 @@ +c +c +c ############################################################ +c ## COPYRIGHT (C) 2022 by Moses KJ Chung & Jay W. Ponder ## +c ## All Rights Reserved ## +c ############################################################ +c +c ################################################################# +c ## ## +c ## subroutine exrepel3 -- exch repulsion energy & analysis ## +c ## ## +c ################################################################# +c +c +c "exrepel3" calculates the exchange repulsion energy and partitions +c the energy among the atoms +c +c literature reference: +c +c M. K. J. Chung and J. W. Ponder, "Beyond isotropic repulsion: +c Classical anisotropic repulsion by inclusion of p orbitals", +c Journal of Chemical Physics, 160, 174118 (2024) +c +c + subroutine exrepel3 + use limits + implicit none +c +c +c choose the method for summing over pairwise interactions +c + if (use_mlist) then + call exrepel3b + else + call exrepel3a + end if + return + end +c +c +c ################################################################## +c ## ## +c ## subroutine exrepel3a -- exch repulsion analysis via loop ## +c ## ## +c ################################################################## +c +c +c "exrepel3a" calculates the exchange repulsion energy and also +c partitions the energy among the atoms using a double loop +c +c + subroutine exrepel3a + use action + use analyz + use atomid + use atoms + use bound + use cell + use couple + use energi + use group + use inform + use inter + use iounit + use molcul + use mutant + use reppot + use shunt + use units + use usage + use xrepel + implicit none + integer i,j,k + integer ii,kk + integer ind1,ind2,ind3 + integer jcell + real*8 e,fgrp,taper + real*8 xi,yi,zi + real*8 xr,yr,zr + real*8 r,r2,r3,r4,r5 + real*8 normi + real*8 zxri,zxrk + real*8 vali,valk + real*8 dmpi,dmpk + real*8 dis,dks + real*8 dmpip,dmpkp + real*8 cis,cks + real*8 cix,ckx + real*8 ciy,cky + real*8 ciz,ckz + real*8 rcix,rckx + real*8 rciy,rcky + real*8 rciz,rckz + real*8 cscs,cxcx,cycy + real*8 czcz,cscz,czcs + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 intS,intS2 + real*8 bi(3) + real*8 bj(3) + real*8 bk(3) + real*8, allocatable :: rscale(:) + logical proceed,usei + logical muti,mutk,mutik + logical header,huge + logical grad + character*6 mode +c +c +c zero out the repulsion energy and partitioning terms +c + ner = 0 + er = 0.0d0 + do i = 1, n + aer(i) = 0.0d0 + end do + if (nxrep .eq. 0) return +c +c check the sign of multipole components at chiral sites +c + call chkpole +c +c determine pseudo orbital coefficients +c + call solvcoeff +c +c rotate the coefficient components into the global frame +c + call rotcoeff +c +c perform dynamic allocation of some local arrays +c + allocate (rscale(n)) +c +c initialize connected atom exclusion coefficients +c + do i = 1, n + rscale(i) = 1.0d0 + end do +c +c set cutoff distances and switching coefficients +c + mode = 'REPULS' + call switch (mode) +c +c set gradient mode to false +c + grad = .false. +c +c calculate the exchange repulsion interaction energy term +c + do ii = 1, nxrep-1 + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kk = ii+1, nxrep + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + if (use_bounds) call image (xr,yr,zr) + r2 = xr * xr + yr * yr + zr * zr + if (r2 .le. off2) then + r = sqrt(r2) + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, 0.0d0, + & r, grad, SS, dSS, SPz, dSPz, PzS, dPzS, + & PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + e = hartree*(zxri*valk+zxrk*vali)*intS2/r*rscale(k) +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + r3 = r2 * r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + e = e * taper + end if +c +c scale the interaction based on its group membership +c + if (use_group) e = e * fgrp +c +c increment the overall exchange repulsion energy component +c + if (e .ne. 0.0d0) then + ner = ner + 1 + er = er + e + aer(i) = aer(i) + 0.5d0*e + aer(k) = aer(k) + 0.5d0*e + end if +c +c increment the total intermolecular energy +c + if (molcule(i) .ne. molcule(k)) then + einter = einter + e + end if + end if + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do +c +c for periodic boundary conditions with large cutoffs +c neighbors must be found by the replicates method +c + if (use_replica) then +c +c calculate interaction energy with other unit cells +c + do ii = 1, nxrep + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kk = ii, nxrep + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + do jcell = 2, ncell + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + call imager (xr,yr,zr,jcell) + r2 = xr*xr + yr* yr + zr*zr + if (r2 .le. off2) then + r = sqrt(r2) + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, + & 0.0d0, r, grad, SS, dSS, SPz, dSPz, PzS, + & dPzS, PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + e = hartree*(zxri*valk+zxrk*vali)*intS2/r* + & rscale(k) +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + r3 = r2 * r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + e = e * taper + end if +c +c scale the interaction based on its group membership +c + if (use_group) e = e * fgrp +c +c increment the overall exchange repulsion energy component; +c interaction of an atom with its own image counts half +c + if (e .ne. 0.0d0) then + ner = ner + 1 + if (i .eq. k) then + er = er + 0.5d0*e + aer(i) = aer(i) + 0.5d0*e + else + er = er + e + aer(i) = aer(i) + 0.5d0*e + aer(k) = aer(k) + 0.5d0*e + end if + end if +c +c increment the total intermolecular energy +c + einter = einter + e + end if + end do + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do + end if +c +c perform deallocation of some local arrays +c + deallocate (rscale) + return + end +c +c +c ################################################################## +c ## ## +c ## subroutine exrepel3b -- exch repulsion analysis via list ## +c ## ## +c ################################################################## +c +c +c "exrepel3b" calculates the exchange repulsion energy and also +c partitions the energy among the atoms using a neighbor list +c +c + subroutine exrepel3b + use action + use analyz + use atomid + use atoms + use bound + use couple + use energi + use group + use inform + use inter + use iounit + use molcul + use mutant + use neigh + use reppot + use shunt + use units + use usage + use xrepel + implicit none + integer i,j,k + integer ii,kk,kkk + integer ind1,ind2,ind3 + real*8 e,fgrp,taper + real*8 xi,yi,zi + real*8 xr,yr,zr + real*8 r,r2,r3,r4,r5 + real*8 normi + real*8 zxri,zxrk + real*8 vali,valk + real*8 dmpi,dmpk + real*8 dis,dks + real*8 dmpip,dmpkp + real*8 cis,cks + real*8 cix,ckx + real*8 ciy,cky + real*8 ciz,ckz + real*8 rcix,rckx + real*8 rciy,rcky + real*8 rciz,rckz + real*8 cscs,cxcx,cycy + real*8 czcz,cscz,czcs + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 intS,intS2 + real*8 bi(3) + real*8 bj(3) + real*8 bk(3) + real*8, allocatable :: rscale(:) + logical proceed,usei + logical muti,mutk,mutik + logical header,huge + logical grad + character*6 mode +c +c +c zero out the repulsion energy and partitioning terms +c + ner = 0 + er = 0.0d0 + do i = 1, n + aer(i) = 0.0d0 + end do + if (nxrep .eq. 0) return +c +c check the sign of multipole components at chiral sites +c + call chkpole +c +c determine pseudo orbital coefficients +c + call solvcoeff +c +c rotate the coefficient components into the global frame +c + call rotcoeff +c +c perform dynamic allocation of some local arrays +c + allocate (rscale(n)) +c +c initialize connected atom exclusion coefficients +c + do i = 1, n + rscale(i) = 1.0d0 + end do +c +c set cutoff distances and switching coefficients +c + mode = 'REPULS' + call switch (mode) +c +c set gradient mode to false +c + grad = .false. +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(nxrep,ixrep,x,y,z,zpxr,dmppxr,rcpxr,n12,i12, +!$OMP& n13,i13,n14,i14,n15,i15,r2scale,r3scale,r4scale,r5scale, +!$OMP& nelst,elst,use,use_group,use_intra,use_bounds,grad, +!$OMP& mut,cut2,off2,c0,c1,c2,c3,c4,c5,molcule,name, +!$OMP& verbose,debug,header,iout) +!$OMP& firstprivate(rscale) +!$OMP& shared (er,ner,aer,einter) +!$OMP DO reduction(+:er,ner,aer,einter) schedule(guided) +c +c calculate the exchange repulsion interaction energy term +c + do ii = 1, nxrep + i = ixrep(ii) + xi = x(i) + yi = y(i) + zi = z(i) + zxri = zpxr(i) + dmpi = dmppxr(i) + vali = zxri + dis = 1.0d0 + dmpip = dis * dmpi + cis = rcpxr(1,i) + cix = rcpxr(2,i) + ciy = rcpxr(3,i) + ciz = rcpxr(4,i) + usei = use(i) + muti = mut(i) +c +c set exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = r2scale + end do + do j = 1, n13(i) + rscale(i13(j,i)) = r3scale + end do + do j = 1, n14(i) + rscale(i14(j,i)) = r4scale + end do + do j = 1, n15(i) + rscale(i15(j,i)) = r5scale + end do +c +c evaluate all sites within the cutoff distance +c + do kkk = 1, nelst(ii) + kk = elst(kkk,ii) + k = ixrep(kk) + mutk = mut(k) + proceed = .true. + if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) + if (.not. use_intra) proceed = .true. + if (proceed) proceed = (usei .or. use(k)) + if (proceed) then + xr = x(k) - xi + yr = y(k) - yi + zr = z(k) - zi + if (use_bounds) call image (xr,yr,zr) + r2 = xr * xr + yr * yr + zr * zr + if (r2 .le. off2) then + r = sqrt(r2) + zxrk = zpxr(k) + dmpk = dmppxr(k) + valk = zxrk + dks = 1.0d0 + dmpkp = dks * dmpk + cks = rcpxr(1,k) + ckx = rcpxr(2,k) + cky = rcpxr(3,k) + ckz = rcpxr(4,k) +c +c choose orthogonal 2-body coordinates / solve rotation matrix +c + bk(1) = xr / r + bk(2) = yr / r + bk(3) = zr / r + ind1 = maxloc(abs(bk), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + bi(ind1) = -bk(ind2) + bi(ind2) = bk(ind1) + bi(ind3) = 0.0d0 + normi = sqrt(bi(1)**2 + bi(2)**2 + bi(3)**2) + bi(1) = bi(1) / normi + bi(2) = bi(2) / normi + bi(3) = bi(3) / normi + bj(1) = bk(2)*bi(3) - bk(3)*bi(2) + bj(2) = bk(3)*bi(1) - bk(1)*bi(3) + bj(3) = bk(1)*bi(2) - bk(2)*bi(1) +c +c rotate p orbital cofficients to 2-body (prolate spheroid) frame +c + rcix = bi(1)*cix + bi(2)*ciy + bi(3)*ciz + rciy = bj(1)*cix + bj(2)*ciy + bj(3)*ciz + rciz = bk(1)*cix + bk(2)*ciy + bk(3)*ciz + rckx = bi(1)*ckx + bi(2)*cky + bi(3)*ckz + rcky = bj(1)*ckx + bj(2)*cky + bj(3)*ckz + rckz = bk(1)*ckx + bk(2)*cky + bk(3)*ckz + cscs = cis * cks + cxcx = rcix * rckx + cycy = rciy * rcky + czcz = rciz * rckz + cscz = cis * rckz + czcs = rciz * cks + call computeOverlap (dmpi, dmpk, dmpip, dmpkp, 0.0d0, + & r, grad, SS, dSS, SPz, dSPz, PzS, dPzS, + & PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + intS = cscs * SS + cxcx * PxPx + cycy * PyPy + & + czcz * PzPz + cscz * SPz + czcs * PzS + intS2 = intS * intS + e = hartree*(zxri*valk+zxrk*vali)*intS2/r*rscale(k) +c +c use energy switching if near the cutoff distance +c + if (r2 .gt. cut2) then + r3 = r2 * r + r4 = r2 * r2 + r5 = r2 * r3 + taper = c5*r5 + c4*r4 + c3*r3 + & + c2*r2 + c1*r + c0 + e = e * taper + end if +c +c scale the interaction based on its group membership +c + if (use_group) e = e * fgrp +c +c increment the overall exchange repulsion energy component +c + if (e .ne. 0.0d0) then + ner = ner + 1 + er = er + e + aer(i) = aer(i) + 0.5d0*e + aer(k) = aer(k) + 0.5d0*e + end if +c +c increment the total intermolecular energy +c + if (molcule(i) .ne. molcule(k)) then + einter = einter + e + end if + end if + end if + end do +c +c reset exclusion coefficients for connected atoms +c + do j = 1, n12(i) + rscale(i12(j,i)) = 1.0d0 + end do + do j = 1, n13(i) + rscale(i13(j,i)) = 1.0d0 + end do + do j = 1, n14(i) + rscale(i14(j,i)) = 1.0d0 + end do + do j = 1, n15(i) + rscale(i15(j,i)) = 1.0d0 + end do + end do +c +c OpenMP directives for the major loop structure +c +!$OMP END DO +!$OMP END PARALLEL +c +c perform deallocation of some local arrays +c + deallocate (rscale) + return + end +c +c +c ################################################################ +c ## ## +c ## subroutine solvcoeff -- solve for orbital coefficients ## +c ## ## +c ################################################################ +c +c +c "solvcoeff" finds the coefficients for the pseudo orbital +c +c + subroutine solvcoeff + use repel + use xrepel + implicit none + integer ii,k + integer ind1,ind2,ind3 + real*8 cr,cs + real*8 pcoeff(3) + real*8 ppole(3) + real*8 p2p1,p3p1 + logical l1,l2,l3 +c +c +c determine pseudo orbital coefficients +c + do ii = 1, nxrep + do k = 1, 3 + pcoeff(k) = 0.0d0 + end do + ppole(1) = xrepole(2,ii) + ppole(2) = xrepole(3,ii) + ppole(3) = xrepole(4,ii) + cr = crpxr(ii) + l1 = (abs(ppole(1)) < 1.0d-10) + l2 = (abs(ppole(2)) < 1.0d-10) + l3 = (abs(ppole(3)) < 1.0d-10) +c +c case for no dipole +c + if (l1.and.l2.and.l3) then + cs = 1.0d0 + ind1 = 1 + ind2 = 2 + ind3 = 3 +c +c case for p orbital coefficients set to 0 +c + else if (cr < 1.0d-10) then + cs = 1.0d0 + ind1 = 1 + ind2 = 2 + ind3 = 3 +c +c determine normalized coefficients +c + else + cs = 1.0d0 / sqrt(1.0d0 + cr) +c +c determine index for largest absolute dipole component +c + ind1 = maxloc(abs(ppole), dim=1) + ind2 = mod(ind1,3) + 1 + ind3 = mod(ind1+1,3) + 1 + p2p1 = ppole(ind2) / ppole(ind1) + p3p1 = ppole(ind3) / ppole(ind1) + pcoeff(ind1) = cs * sqrt(cr / (1.0d0 + p2p1**2 + p3p1**2)) + if (ppole(ind1) < 0.0d0) then + pcoeff(ind1) = -pcoeff(ind1) + end if + pcoeff(ind2) = pcoeff(ind1) * p2p1 + pcoeff(ind3) = pcoeff(ind1) * p3p1 + end if + cpxr(1,ii) = cs + cpxr(ind1+1,ii) = pcoeff(ind1) + cpxr(ind2+1,ii) = pcoeff(ind2) + cpxr(ind3+1,ii) = pcoeff(ind3) + end do + return + end +c +c +c ############################################################ +c ## ## +c ## subroutine rotcoeff -- rotate orbital coefficients ## +c ## ## +c ############################################################ +c +c +c "rotcoeff" rotates the coefficients for the pseudo orbital +c to global coordinates +c +c + subroutine rotcoeff + use mpole + use repel + use xrepel + implicit none + integer isite + integer i,j + real*8 a(3,3) + real*8 cpole(4) + logical planar + character*8 axetyp +c +c +c rotate pseudo orbital coefficients +c + do isite = 1, nxrep +c +c determine rotation matrix +c + call rotmat (isite,a,planar) +c +c copy local frame coefficients +c + do i = 1, 4 + cpole(i) = cpxr(i,isite) + end do + if (planar) then + axetyp = polaxe(isite) + if (axetyp .eq. 'Z-Bisect') then + cpole(2) = 0.0d0 + else if (axetyp .eq. '3-Fold') then + do i = 2, 4 + cpole(i) = 0.0d0 + end do + end if + end if +c +c s orbital coefficients have the same value in any coordinate frame +c + rcpxr(1,isite) = cpole(1) +c +c rotate the p orbital coefficients to the global coordinate frame +c + do i = 2, 4 + rcpxr(i,isite) = 0.0d0 + do j = 2, 4 + rcpxr(i,isite) = rcpxr(i,isite) + cpole(j)*a(i-1,j-1) + end do + end do + end do + return + end +c +c +c ################################################################# +c ## ## +c ## subroutine computeOverlap -- computes overlap integrals ## +c ## ## +c ################################################################# +c +c +c "computeOverlap" computes overlap integrals +c +c + subroutine computeOverlap (a, b, ap, bp, z1, z2, grad, SS, dSS, + & SPz, dSPz, PzS, dPzS, PxPx, dPxPx, PyPy, dPyPy, PzPz, dPzPz) + implicit none + real*8 a,b,ap,bp + real*8 z1,z2 + real*8 SS,SPz,PzS + real*8 PxPx,PyPy,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPyPy,dPzPz + real*8 eps, diffa, diffb + logical grad +c +c + eps = 1.0d-10 + diffa = abs(a - ap) + diffb = abs(b - bp) + if (diffa.lt.eps .and. diffb.lt.eps) then + call overlapAll(a, b, z1, z2, grad, SS, dSS, SPz, dSPz, + & PzS, dPzS, PxPx, dPxPx, PzPz, dPzPz) + else + call overlapSS(a, b, z1, z2, grad, SS, dSS) + call overlapSPz(a, bp, z1, z2, grad, SPz, dSPz) + call overlapPzS(ap, b, z1, z2, grad, PzS, dPzS) + call overlapPxPx(ap, bp, z1, z2, grad, PxPx, dPxPx) + call overlapPzPz(ap, bp, z1, z2, grad, PzPz, dPzPz) + end if + PyPy = PxPx + dPyPy = dPxPx + return + end +c +c +c ######################################################## +c ## ## +c ## subroutine overlapAll -- all overlap integrals ## +c ## ## +c ######################################################## +c +c +c "overlapAll" computes all the overlap integrals +c +c + subroutine overlapAll (a, b, z1, z2, grad, SS, dSS, SPz, dSPz, + & PzS, dPzS, PxPx, dPxPx, PzPz, dPzPz) + implicit none + real*8 SS,SPz,PzS + real*8 PxPx,PzPz + real*8 dSS,dSPz,dPzS + real*8 dPxPx,dPzPz + real*8 a,b,z1,z2 + real*8 diff,eps,zr,r + real*8 rho,rho2,rho3,rho4 + real*8 exp1 + real*8 alpha,tau,tau2 + real*8 rhoA,rhoB + real*8 a2,b2,kappa + real*8 pre,pre1,pre2,pre3 + real*8 term1,term2 + real*8 rhoA2,rhoA3,rhoA4,rhoA5 + real*8 rhoB2,rhoB3,rhoB4,rhoB5 + real*8 kappam,kappam2 + real*8 kappap,kappap2 + real*8 taurho,taurho2,taurho3 + real*8 expA,expB + logical grad +c +c + diff = abs(a - b) + eps = 0.001d0 + zr = z2 - z1 + r = abs(zr) + if (diff .lt. eps) then + rho = a * r + rho2 = rho * rho + rho3 = rho2 * rho + rho4 = rho3 * rho + exp1 = exp(-rho) + SS = (1.0d0 + rho + rho2 / 3.0d0) * exp1 + SPz = -0.5d0 * rho * (1.0d0 + rho + rho2 / 3.0d0) * exp1 + PzS = -SPz + PxPx = (1.0d0 + rho + 2.0d0/5.0d0 * rho2 + rho3/15.0d0) * exp1 + PzPz = -(-1.0d0 - rho - rho2 / 5.0d0 + 2.0d0/15.0d0 * rho3 + & + rho4 / 15.0d0) * exp1 + if (grad) then + dSS = -1.0d0/3.0d0 * a * rho * (1.0d0 + rho) * exp1 + dSPz = -0.5d0 * a * (1.0d0 + rho - rho3 / 3.0d0) * exp1 + dPzS = -dSPz + dPxPx = -0.2d0 * a * rho * (1.0d0 + rho + rho2 / 3.0d0)*exp1 + dPzPz = -0.6d0 * a * rho * (1.0d0 + rho + 2.0d0/9.0d0 * rho2 + & - rho3 / 9.0d0) * exp1 + end if + else + alpha = 1.0d0 / 2.0d0 * (a + b) + tau = (a - b) / (a + b) + tau2 = tau * tau + rho = alpha * r + rhoA = a * r + rhoB = b * r + a2 = a * a + b2 = b * b + kappa = (a2 + b2) / (a2 - b2) + rho2 = rho * rho + rho3 = rho2 * rho + rhoA2 = rhoA * rhoA + rhoA3 = rhoA2 * rhoA + rhoA4 = rhoA3 * rhoA + rhoB2 = rhoB * rhoB + rhoB3 = rhoB2 * rhoB + rhoB4 = rhoB3 * rhoB + kappam = 1.0d0 - kappa + kappap = 1.0d0 + kappa + kappam2 = kappam**2 + kappap2 = kappap**2 + taurho = tau * rho + taurho2 = taurho * rho + taurho3 = taurho2 * rho + expA = exp(-rhoA) + expB = exp(-rhoB) + pre1 = sqrt(1.0d0 - tau2) + pre2 = sqrt((1.0d0 + tau) / (1.0d0 - tau)) + pre3 = sqrt((1.0d0 - tau) / (1.0d0 + tau)) + pre = pre1 / taurho + term1 =-kappam * (2.0d0 * kappap + rhoA) * expA + term2 = kappap * (2.0d0 * kappam + rhoB) * expB + SS = pre * (term1 + term2) + pre = pre2 / taurho2 + term1 =-kappam2 * (6.0d0 * kappap * (1.0d0 + rhoA) + & + 2.0d0 * rhoA2) * expA + term2 = kappap * (6.0d0 * kappam2 * (1.0d0 + rhoB) + & + 4.0d0 * kappam * rhoB2 + rhoB3) * expB + SPz = -pre * (term1 + term2) + pre = -pre3 / taurho2 + term1 =-kappap2 * (6.0d0 * kappam * (1.0d0 + rhoB) + & + 2.0d0 * rhoB2) * expB + term2 = kappam * (6.0d0 * kappap2 * (1.0d0 + rhoA) + & + 4.0d0 * kappap * rhoA2 + rhoA3) * expA + PzS = pre * (term1 + term2) + pre = 1.0d0 / (pre1 * taurho3) + term1 =-kappam2 * (24.0d0 * kappap2 * (1.0d0 + rhoA) + & + 12.0d0 * kappap * rhoA2 + 2.0d0 * rhoA3) * expA + term2 = kappap2 * (24.0d0 * kappam2 * (1.0d0 + rhoB) + & + 12.0d0 * kappam * rhoB2 + 2.0d0 * rhoB3) * expB + PxPx = pre * (term1 + term2) + term1 =-kappam2 * (48.0d0 * kappap2 * (1.0d0 + rhoA + & + 0.5d0 * rhoA2) + 2.0d0 * (5.0d0 + 6.0d0 * kappa) * rhoA3 + & + 2.0d0 * rhoA4) * expA + term2 = kappap2 * (48.0d0 * kappam2 * (1.0d0 + rhoB + & + 0.5d0 * rhoB2) + 2.0d0 * (5.0d0 - 6.0d0 * kappa) * rhoB3 + & + 2.0d0 * rhoB4) * expB + PzPz = -pre * (term1 + term2) + if (grad) then + rhoA5 = rhoA4 * rhoA + rhoB5 = rhoB4 * rhoB + pre = pre1 / taurho + term1 = kappam * (2.0d0 * kappap * (1.0d0 + rhoA) + rhoA2) + & * expA + term2 =-kappap * (2.0d0 * kappam * (1.0d0 + rhoB) + rhoB2) + & * expB + dSS = pre / r * (term1 + term2) + pre = pre2 / taurho2 + term1 = 2.0d0 * kappam2 * (6.0d0 * kappap * (1.0d0 + rhoA + & + 0.5d0 * rhoA2) + rhoA3) * expA + term2 = kappap * (-12.0d0 * kappam2 * (1.0d0 + rhoB + 0.5d0 + & * rhoB2) + (1.0d0 - 4.0d0 * kappam) * rhoB3 - rhoB4) + & * expB + dSPz = -pre / r * (term1 + term2) + pre = -pre3 / taurho2 + term1 = 2.0d0 * kappap2 * (6.0d0 * kappam * (1.0d0 + rhoB + & + 0.5d0 * rhoB2) + rhoB3) * expB + term2 = kappam * (-12.0d0 * kappap2 * (1.0d0 + rhoA + 0.5d0 + & * rhoA2) + (1.0d0 - 4.0d0 * kappap) * rhoA3 - rhoA4) + & * expA + dPzS = pre / r * (term1 + term2) + pre = 1.0d0 / (pre1 * taurho3) + term1 = 2.0d0 * kappam2 * (36.0d0 * kappap2 * (1.0d0 + rhoA) + & + 6.0d0 * kappap * (1.0d0 + 2.0d0 * kappap) * rhoA2 + & + 6.0d0 * kappap * rhoA3 + rhoA4) * expA + term2 =-2.0d0 * kappap2 * (36.0d0 * kappam2 * (1.0d0 + rhoB) + & + 6.0d0 * kappam * (1.0d0 + 2.0d0 * kappam) * rhoB2 + & + 6.0d0 * kappam * rhoB3 + rhoB4) * expB + dPxPx = pre / r * (term1 + term2) + term1 = kappam2 * (72.0d0 * kappap2 * (1.0d0 + rhoA + & + 0.5d0 * rhoA2 + rhoA3 / 6.0d0) + 2.0d0 * (2.0d0 + & + 3.0d0 * kappa) * rhoA4 + rhoA5) * expA + term2 =-kappap2 * (72.0d0 * kappam2 * (1.0d0 + rhoB + & + 0.5d0 * rhoB2 + rhoB3 / 6.0d0) + 2.0d0 * (2.0d0 + & - 3.0d0 * kappa) * rhoB4 + rhoB5) * expB + dPzPz = -2.0d0 * pre / r * (term1 + term2) + end if + end if + return + end +c +c +c ##################################################### +c ## ## +c ## subroutine overlapSS -- SS overlap integral ## +c ## ## +c ##################################################### +c +c +c "overlapSS" computes the overlap integral +c +c + subroutine overlapSS (a, b, z1, z2, grad, SS, dSS) + implicit none + real*8 over,dover + real*8 SS,dSS + real*8 a,b,z1,z2 + real*8 diff,eps,zr,r + real*8 rho,rho2 + real*8 exp1 + real*8 alpha,tau,rhoA,rhoB + real*8 a2,b2,kappa + real*8 pre,term1,term2 + real*8 rhoA2,rhoB2 + real*8 kappam,kappap + real*8 expA,expB + logical grad +c +c + diff = abs(a - b) + eps = 0.001d0 + zr = z2 - z1 + r = abs(zr) + if (diff .lt. eps) then + rho = a * r + rho2 = rho * rho + exp1 = exp(-rho) + over = (1.0d0 + rho + rho2 / 3.0d0) * exp1 + if (grad) then + dover = -1.0d0/3.0d0 * a * rho * (1.0d0 + rho) * exp1 + end if + else + alpha = 1.0d0 / 2.0d0 * (a + b) + tau = (a - b) / (a + b) + rho = alpha * r + rhoA = a * r + rhoB = b * r + a2 = a * a + b2 = b * b + kappa = (a2 + b2) / (a2 - b2) + kappam = 1.0d0 - kappa + kappap = 1.0d0 + kappa + expA = exp(-rhoA) + expB = exp(-rhoB) + pre = sqrt(1.0d0 - tau**2) / (tau * rho) + term1 =-kappam * (2.0d0 * kappap + rhoA) * expA + term2 = kappap * (2.0d0 * kappam + rhoB) * expB + over = pre * (term1 + term2) + if (grad) then + rhoA2 = rhoA * rhoA + rhoB2 = rhoB * rhoB + term1 = kappam * (2.0d0 * kappap * (1.0d0 + rhoA) + rhoA2) + & * expA + term2 =-kappap * (2.0d0 * kappam * (1.0d0 + rhoB) + rhoB2) + & * expB + dover = pre / r * (term1 + term2) + end if + end if + SS = over + dSS = dover + return + end +c +c +c ####################################################### +c ## ## +c ## subroutine overlapSPz -- SPz overlap integral ## +c ## ## +c ####################################################### +c +c +c "overlapSPz" computes the overlap integral +c +c + subroutine overlapSPz (a, b, z1, z2, grad, SPz, dSPz) + implicit none + real*8 over,dover + real*8 SPz,dSPz + real*8 a,b,z1,z2 + real*8 diff,eps,zr,r + real*8 rho,rho2,rho3 + real*8 exp1 + real*8 alpha,tau,rhoA,rhoB + real*8 a2,b2,kappa + real*8 pre,term1,term2 + real*8 rhoA2,rhoA3 + real*8 rhoB2,rhoB3,rhoB4 + real*8 kappam,kappam2 + real*8 kappap,kappap2 + real*8 expA,expB + logical grad +c +c + diff = abs(a - b) + eps = 0.001d0 + zr = z2 - z1 + r = abs(zr) + if (diff .lt. eps) then + rho = a * r + rho2 = rho * rho + exp1 = exp(-rho) + over = 0.5d0 * rho * (1.0d0 + rho + rho2 / 3.0d0) * exp1 + if (grad) then + rho3 = rho2 * rho + dover = 0.5d0 * a * (1.0d0 + rho - rho3 / 3.0d0) * exp1 + end if + else + alpha = 1.0d0 / 2.0d0 * (a + b) + tau = (a - b) / (a + b) + rho = alpha * r + rhoA = a * r + rhoB = b * r + a2 = a * a + b2 = b * b + kappa = (a2 + b2) / (a2 - b2) + rho2 = rho**2 + rhoA2 = rhoA * rhoA + rhoB2 = rhoB * rhoB + rhoB3 = rhoB2 * rhoB + kappam = 1.0d0 - kappa + kappap = 1.0d0 + kappa + kappam2 = kappam**2 + kappap2 = kappap**2 + expA = exp(-rhoA) + expB = exp(-rhoB) + pre = sqrt((1.0d0 + tau) / (1.0d0 - tau)) / (tau * rho2) + term1 =-kappam2 * (6.0d0 * kappap * (1.0d0 + rhoA) + & + 2.0d0 * rhoA2) * expA + term2 = kappap * (6.0d0 * kappam2 * (1.0d0 + rhoB) + & + 4.0d0 * kappam * rhoB2 + rhoB3) * expB + over = pre * (term1 + term2) + if (grad) then + rhoA3 = rhoA2 * rhoA + rhoB4 = rhoB3 * rhoB + term1 = 2.0d0 * kappam2 * (6.0d0 * kappap * (1.0d0 + rhoA + & + 0.5d0 * rhoA2) + rhoA3) * expA + term2 = kappap * (-12.0d0 * kappam2 * (1.0d0 + rhoB + 0.5d0 + & * rhoB2) + (1.0d0 - 4.0d0 * kappam) * rhoB3 - rhoB4) + & * expB + dover = pre / r * (term1 + term2) + end if + end if + SPz = -over + dSPz = -dover + if (z1 > z2) then + SPz = -SPz + dSPz = -dSPz + end if + return + end +c +c +c ####################################################### +c ## ## +c ## subroutine overlapPzS -- PzS overlap integral ## +c ## ## +c ####################################################### +c +c +c "overlapPzS" computes the overlap integral +c +c + subroutine overlapPzS (a, b, z1, z2, grad, PzS, dPzS) + implicit none + real*8 PzS,dPzS + real*8 a,b,z1,z2 + logical grad +c +c + call overlapSPz(b, a, z2, z1, grad, PzS, dPzS) + return + end +c +c +c ######################################################### +c ## ## +c ## subroutine overlapPxPx -- PxPx overlap integral ## +c ## ## +c ######################################################### +c +c +c "overlapPxPx" computes the overlap integral +c +c + subroutine overlapPxPx (a, b, z1, z2, grad, PxPx, dPxPx) + implicit none + real*8 over,dover + real*8 PxPx,dPxPx + real*8 a,b,z1,z2 + real*8 diff,eps,zr,r + real*8 rho,rho2,rho3 + real*8 exp1 + real*8 alpha,tau,rhoA,rhoB + real*8 a2,b2,kappa + real*8 pre,term1,term2 + real*8 rhoA2,rhoA3,rhoA4 + real*8 rhoB2,rhoB3,rhoB4 + real*8 kappam,kappam2 + real*8 kappap,kappap2 + real*8 expA,expB + logical grad +c +c + diff = abs(a - b) + eps = 0.001d0 + zr = z2 - z1 + r = abs(zr) + if (diff .lt. eps) then + rho = a * r + rho2 = rho * rho + rho3 = rho2 * rho + exp1 = exp(-rho) + over = (1.0d0 + rho + 2.0d0/5.0d0 * rho2 + rho3/15.0d0) * exp1 + if (grad) then + dover = -0.2d0 * a * rho * (1.0d0 + rho + rho2 / 3.0d0) + & * exp1 + end if + else + alpha = 1.0d0 / 2.0d0 * (a + b) + tau = (a - b) / (a + b) + rho = alpha * r + rhoA = a * r + rhoB = b * r + a2 = a * a + b2 = b * b + kappa = (a2 + b2) / (a2 - b2) + rho3 = rho**3 + rhoA2 = rhoA * rhoA + rhoA3 = rhoA2 * rhoA + rhoB2 = rhoB * rhoB + rhoB3 = rhoB2 * rhoB + kappam = 1.0d0 - kappa + kappap = 1.0d0 + kappa + kappam2 = kappam**2 + kappap2 = kappap**2 + expA = exp(-rhoA) + expB = exp(-rhoB) + pre = 1.0d0 / (sqrt(1.0d0 - tau**2) * tau * rho3) + term1 =-kappam2 * (24.0d0 * kappap2 * (1.0d0 + rhoA) + & + 12.0d0 * kappap * rhoA2 + 2.0d0 * rhoA3) * expA + term2 = kappap2 * (24.0d0 * kappam2 * (1.0d0 + rhoB) + & + 12.0d0 * kappam * rhoB2 + 2.0d0 * rhoB3) * expB + over = pre * (term1 + term2) + if (grad) then + rhoA4 = rhoA3 * rhoA + rhoB4 = rhoB3 * rhoB + term1 = 2.0d0 * kappam2 * (36.0d0 * kappap2 * (1.0d0 + rhoA) + & + 6.0d0 * kappap * (1.0d0 + 2.0d0 * kappap) * rhoA2 + & + 6.0d0 * kappap * rhoA3 + rhoA4) * expA + term2 =-2.0d0 * kappap2 * (36.0d0 * kappam2 * (1.0d0 + rhoB) + & + 6.0d0 * kappam * (1.0d0 + 2.0d0 * kappam) * rhoB2 + & + 6.0d0 * kappam * rhoB3 + rhoB4) * expB + dover = pre / r * (term1 + term2) + end if + end if + PxPx = over + dPxPx = dover + return + end +c +c +c ######################################################### +c ## ## +c ## subroutine overlapPzPz -- PzPz overlap integral ## +c ## ## +c ######################################################### +c +c +c "overlapPzPz" computes the overlap integral +c +c + subroutine overlapPzPz (a, b, z1, z2, grad, PzPz, dPzPz) + implicit none + real*8 over,dover + real*8 PzPz,dPzPz + real*8 a,b,z1,z2 + real*8 diff,eps,zr,r + real*8 rho,rho2,rho3,rho4 + real*8 exp1 + real*8 alpha,tau,rhoA,rhoB + real*8 a2,b2,kappa + real*8 pre,term1,term2 + real*8 rhoA2,rhoA3,rhoA4,rhoA5 + real*8 rhoB2,rhoB3,rhoB4,rhoB5 + real*8 kappam,kappam2 + real*8 kappap,kappap2 + real*8 expA,expB + logical grad +c +c + diff = abs(a - b) + eps = 0.001d0 + zr = z2 - z1 + r = abs(zr) + if (diff .lt. eps) then + rho = a * r + rho2 = rho * rho + rho3 = rho2 * rho + rho4 = rho3 * rho + exp1 = exp(-rho) + over = (-1.0d0 - rho - rho2 / 5.0d0 + 2.0d0/15.0d0 * rho3 + & + rho4 / 15.0d0) * exp1 + if (grad) then + dover = 0.6d0 * a * rho * (1.0d0 + rho + 2.0d0/9.0d0 * rho2 + & - rho3 / 9.0d0) * exp1 + end if + else + alpha = 1.0d0 / 2.0d0 * (a + b) + tau = (a - b) / (a + b) + rho = alpha * r + rhoA = a * r + rhoB = b * r + a2 = a * a + b2 = b * b + kappa = (a2 + b2) / (a2 - b2) + rho3 = rho**3 + rhoA2 = rhoA * rhoA + rhoA3 = rhoA2 * rhoA + rhoA4 = rhoA3 * rhoA + rhoB2 = rhoB * rhoB + rhoB3 = rhoB2 * rhoB + rhoB4 = rhoB3 * rhoB + kappam = 1.0d0 - kappa + kappap = 1.0d0 + kappa + kappam2 = kappam**2 + kappap2 = kappap**2 + expA = exp(-rhoA) + expB = exp(-rhoB) + pre = 1.0d0 / (sqrt(1.0d0 - tau**2) * tau * rho3) + term1 =-kappam2 * (48.0d0 * kappap2 * (1.0d0 + rhoA + & + 0.5d0 * rhoA2) + 2.0d0 * (5.0d0 + 6.0d0 * kappa) * rhoA3 + & + 2.0d0 * rhoA4) * expA + term2 = kappap2 * (48.0d0 * kappam2 * (1.0d0 + rhoB + & + 0.5d0 * rhoB2) + 2.0d0 * (5.0d0 - 6.0d0 * kappa) * rhoB3 + & + 2.0d0 * rhoB4) * expB + over = pre * (term1 + term2) + if (grad) then + rhoA5 = rhoA4 * rhoA + rhoB5 = rhoB4 * rhoB + term1 = kappam2 * (72.0d0 * kappap2 * (1.0d0 + rhoA + & + 0.5d0 * rhoA2 + rhoA3 / 6.0d0) + 2.0d0 * (2.0d0 + & + 3.0d0 * kappa) * rhoA4 + rhoA5) * expA + term2 =-kappap2 * (72.0d0 * kappam2 * (1.0d0 + rhoB + & + 0.5d0 * rhoB2 + rhoB3 / 6.0d0) + 2.0d0 * (2.0d0 + & - 3.0d0 * kappa) * rhoB4 + rhoB5) * expB + dover = 2.0d0 * pre / r * (term1 + term2) + end if + end if + PzPz = -over + dPzPz = -dover + return + end diff --git a/source/field.f b/source/field.f index de300105..e58220b2 100644 --- a/source/field.f +++ b/source/field.f @@ -50,6 +50,7 @@ subroutine field use_tortor = .true. use_vdw = .true. use_repel = .true. + use_xrepel = .true. use_disp = .true. use_charge = .true. use_chgdpl = .true. diff --git a/source/final.f b/source/final.f index 1f02e0dd..5f40746b 100644 --- a/source/final.f +++ b/source/final.f @@ -83,6 +83,7 @@ subroutine final use kurybr use kvdwpr use kvdws + use kxrepl use light use limits use merck @@ -138,6 +139,7 @@ subroutine final use vdw use vibs use warp + use xrepel implicit none c c @@ -737,6 +739,12 @@ subroutine final if (allocated(eps4)) deallocate (eps4) if (allocated(reduct)) deallocate (reduct) c +c deallocation of global arrays from module kxrepl +c + if (allocated(pxrz)) deallocate (pxrz) + if (allocated(pxrdmp)) deallocate (pxrdmp) + if (allocated(pxrcr)) deallocate (pxrcr) +c c deallocation of global arrays from module light c if (allocated(kbx)) deallocate (kbx) @@ -1055,6 +1063,17 @@ subroutine final if (allocated(repole)) deallocate (repole) if (allocated(rrepole)) deallocate (rrepole) c +c deallocation of global arrays from module xrepel +c + if (allocated(ixrep)) deallocate (ixrep) + if (allocated(xreplist)) deallocate (xreplist) + if (allocated(zpxr)) deallocate (zpxr) + if (allocated(dmppxr)) deallocate (dmppxr) + if (allocated(crpxr)) deallocate (crpxr) + if (allocated(cpxr)) deallocate (cpxr) + if (allocated(rcpxr)) deallocate (rcpxr) + if (allocated(xrepole)) deallocate (xrepole) +c c deallocation of global arrays from module restrn c if (allocated(ipfix)) deallocate (ipfix) diff --git a/source/gradient.f b/source/gradient.f index 3c61c40b..4c64db19 100644 --- a/source/gradient.f +++ b/source/gradient.f @@ -240,7 +240,11 @@ subroutine gradient (energy,derivs) if (vdwtyp .eq. 'BUFFERED-14-7') call ehal1 if (vdwtyp .eq. 'GAUSSIAN') call egauss1 end if - if (use_repel) call erepel1 + if (use_repel) then + call erepel1 + else if (use_xrepel) then + call exrepel1 + end if if (use_disp) call edisp1 c c call any miscellaneous energy and gradient routines diff --git a/source/initprm.f b/source/initprm.f index 0634201b..5b8720ee 100644 --- a/source/initprm.f +++ b/source/initprm.f @@ -57,6 +57,7 @@ subroutine initprm use kurybr use kvdws use kvdwpr + use kxrepl use math use mplpot use polpot @@ -70,6 +71,7 @@ subroutine initprm use units use uprior use vdwpot + use xrepel implicit none integer i,j character*3 blank3 @@ -220,6 +222,9 @@ subroutine initprm if (.not. allocated(prsiz)) allocate (prsiz(maxclass)) if (.not. allocated(prdmp)) allocate (prdmp(maxclass)) if (.not. allocated(prele)) allocate (prele(maxclass)) + if (.not. allocated(pxrz)) allocate (pxrz(maxclass)) + if (.not. allocated(pxrdmp)) allocate (pxrdmp(maxclass)) + if (.not. allocated(pxrcr)) allocate (pxrcr(maxclass)) if (.not. allocated(dspsix)) allocate (dspsix(maxclass)) if (.not. allocated(dspdmp)) allocate (dspdmp(maxclass)) if (.not. allocated(chg)) allocate (chg(maxtyp)) @@ -278,6 +283,9 @@ subroutine initprm prsiz(i) = 0.0d0 prdmp(i) = 0.0d0 prele(i) = 0.0d0 + pxrz(i) = 0.0d0 + pxrdmp(i) = 0.0d0 + pxrcr(i) = 0.0d0 dspsix(i) = 0.0d0 dspdmp(i) = 0.0d0 cpele(i) = 0.0d0 diff --git a/source/kmpole.f b/source/kmpole.f index 46309e66..5778a2ea 100644 --- a/source/kmpole.f +++ b/source/kmpole.f @@ -685,8 +685,8 @@ subroutine kmpole c c remove zero or undefined electrostatic sites from the list c - if ((use_mpole .or. use_repel .or. use_solv) .and. - & .not.use_polar .and. .not.use_chgtrn) then + if ((use_mpole .or. use_repel .or. use_xrepel .or. use_solv) + & .and. .not.use_polar .and. .not.use_chgtrn) then npole = 0 ncp = 0 do i = 1, n diff --git a/source/kpolar.f b/source/kpolar.f index 2ce8df4e..b9d7f884 100644 --- a/source/kpolar.f +++ b/source/kpolar.f @@ -532,7 +532,7 @@ subroutine kpolar c c remove zero or undefined electrostatic sites from the list c - if ((use_polar .or. use_repel .or. use_solv) + if ((use_polar .or. use_repel .or. use_xrepel .or. use_solv) & .and. .not.use_chgtrn) then npole = 0 ncp = 0 diff --git a/source/krepel.f b/source/krepel.f index 7f148759..c5e8f139 100644 --- a/source/krepel.f +++ b/source/krepel.f @@ -138,7 +138,7 @@ subroutine krepel write (iout,50) 50 format (/,' Additional Pauli Repulsion Values', & ' for Specific Atoms :', - & //,8x,'Atom',17x,'Size',12x,'Damp', + & //,8x,'Atom',18x,'Size',11x,'Damp', & 8x,'Valence'/) end if if (.not. silent) then diff --git a/source/kxrepel.f b/source/kxrepel.f new file mode 100644 index 00000000..55c8e18e --- /dev/null +++ b/source/kxrepel.f @@ -0,0 +1,183 @@ +c +c +c ############################################################ +c ## COPYRIGHT (C) 2022 by Moses KJ Chung & Jay W. Ponder ## +c ## All Rights Reserved ## +c ############################################################ +c +c ############################################################## +c ## ## +c ## subroutine kxrepel -- exch repulsion term assignment ## +c ## ## +c ############################################################## +c +c +c "kxrepel" assigns the nuclear charge parameter and exponential +c parameter for exchange repulsion interaction and processes any +c new or changed values for these parameters +c +c + subroutine kxrepel + use atomid + use atoms + use inform + use iounit + use kxrepl + use keys + use mpole + use potent + use xrepel + use sizes + implicit none + integer i,j,k + integer iboys + integer freeunit + integer ia,ic,next + real*8 zpr,apr,cpr + logical header + character*20 keyword + character*240 record + character*240 string +c +c +c process keywords containing exch repulsion parameters +c + header = .true. + do i = 1, nkey + next = 1 + record = keyline(i) + call gettext (record,keyword,next) + call upcase (keyword) + if (keyword(1:11) .eq. 'XREPULSION ') then + k = 0 + zpr = 0.0d0 + apr = 0.0d0 + cpr = 0.0d0 + call getnumb (record,k,next) + string = record(next:240) + read (string,*,err=10,end=10) zpr,apr,cpr + 10 continue + if (k .gt. 0) then + if (header .and. .not.silent) then + header = .false. + write (iout,20) + 20 format (/,' Additional Exchange Repulsion', + & ' Parameters :', + & //,5x,'Atom Class',13x,'Charge',11x,'Damp', + & 7x,'PS-Ratio'/) + end if + if (k .le. maxclass) then + pxrz(k) = zpr + pxrdmp(k) = apr + pxrcr(k) = cpr + if (.not. silent) then + write (iout,30) k,zpr,apr,cpr + 30 format (6x,i6,7x,3f15.4) + end if + else + write (iout,40) + 40 format (/,' KXREPEL -- Too many Exchange Repulsion', + & ' Parameters') + abort = .true. + end if + end if + end if + end do +c +c perform dynamic allocation of some global arrays +c + if (allocated(ixrep)) deallocate (ixrep) + if (allocated(xreplist)) deallocate (xreplist) + if (allocated(zpxr)) deallocate (zpxr) + if (allocated(dmppxr)) deallocate (dmppxr) + if (allocated(crpxr)) deallocate (crpxr) + if (allocated(cpxr)) deallocate (cpxr) + if (allocated(rcpxr)) deallocate (rcpxr) + if (allocated(xrepole)) deallocate (xrepole) + allocate (ixrep(n)) + allocate (xreplist(n)) + allocate (zpxr(n)) + allocate (dmppxr(n)) + allocate (crpxr(n)) + allocate (cpxr(4,n)) + allocate (rcpxr(4,n)) + allocate (xrepole(maxpole,n)) +c +c assign the core, alpha, and coefficient ratio parameters +c + do i = 1, n + ixrep(i) = 0 + xreplist(i) = 0 + zpxr(i) = 0.0d0 + dmppxr(i) = 0.0d0 + crpxr(i) = 0.0d0 + ic = class(i) + if (ic .ne. 0) then + zpxr(i) = pxrz(ic) + dmppxr(i) = pxrdmp(ic) + crpxr(i) = pxrcr(ic) + end if + end do +c +c process keywords containing atom specific exchange repulsion +c + header = .true. + do i = 1, nkey + next = 1 + record = keyline(i) + call gettext (record,keyword,next) + call upcase (keyword) + if (keyword(1:11) .eq. 'XREPULSION ') then + ia = 0 + zpr = 0.0d0 + apr = 0.0d0 + cpr = 0.0d0 + string = record(next:240) + read (string,*,err=70,end=70) ia,zpr,apr,cpr + if (ia.lt.0 .and. ia.ge.-n) then + ia = -ia + if (header .and. .not.silent) then + header = .false. + write (iout,50) + 50 format (/,' Additional Exchange Repulsion Values', + & ' for Specific Atoms :', + & //,8x,'Atom',16x,'Charge',11x,'Damp', + & 7x,'PS-Ratio'/) + end if + if (.not. silent) then + write (iout,60) ia,zpr,apr,cpr + 60 format (6x,i6,7x,3f15.4) + end if + zpxr(ia) = zpr + dmppxr(ia) = apr + crpxr(ia) = cpr + end if + 70 continue + end if + end do +c +c condense repulsion sites to the list of multipole sites +c + nxrep = 0 + if (use_xrepel) then + do i = 1, n + if (zpxr(i) .ne. 0) then + nxrep = nxrep + 1 + ixrep(nxrep) = i + xreplist(i) = nxrep + do j = 1, maxpole + xrepole(j,i) = pole(j,i) + end do + end if + end do + end if +c +c test multipoles at chiral sites and invert if necessary +c + call chkpole +c +c turn off the exchange repulsion potential if not used +c + if (nxrep .eq. 0) use_xrepel = .false. + return + end diff --git a/source/kxrepl.f b/source/kxrepl.f new file mode 100644 index 00000000..661c6b74 --- /dev/null +++ b/source/kxrepl.f @@ -0,0 +1,26 @@ +c +c +c ############################################################ +c ## COPYRIGHT (C) 2022 by Moses KJ Chung & Jay W. Ponder ## +c ## All Rights Reserved ## +c ############################################################ +c +c ############################################################### +c ## ## +c ## module kxrepl -- exch repulsion forcefield parameters ## +c ## ## +c ############################################################### +c + +c pxrz nuclear charge parameter value for each atom class +c pxrdmp exch repulsion alpha parameter for each atom class +c pxrcr ratio of p/s orbital cofficients for each atom class +c +c + module kxrepl + implicit none + real*8, allocatable :: pxrz(:) + real*8, allocatable :: pxrdmp(:) + real*8, allocatable :: pxrcr(:) + save + end diff --git a/source/mechanic.f b/source/mechanic.f index 16f46b27..838b12fc 100644 --- a/source/mechanic.f +++ b/source/mechanic.f @@ -107,6 +107,7 @@ subroutine mechanic c call kvdw call krepel + call kxrepel call kdisp c c assign solvation, metal, pisystem and restraint parameters diff --git a/source/nblist.f b/source/nblist.f index cd31c9df..5a51c063 100644 --- a/source/nblist.f +++ b/source/nblist.f @@ -42,7 +42,7 @@ subroutine nblist if (use_vdw .and. use_vlist) call vlist if (use_disp .and. use_dlist) call dlist if ((use_charge.or.use_solv) .and. use_clist) call clist - if ((use_repel.or.use_mpole.or.use_polar + if ((use_repel.or.use_xrepel.or.use_mpole.or.use_polar & .or.use_chgtrn.or.use_solv) .and. use_mlist) call mlist if (use_polar .and. use_ulist) call ulist return diff --git a/source/potent.f b/source/potent.f index 7c2bebe2..8098de3e 100644 --- a/source/potent.f +++ b/source/potent.f @@ -28,6 +28,7 @@ c use_tortor logical flag governing use of torsion-torsion term c use_vdw logical flag governing use of van der Waals potential c use_repel logical flag governing use of Pauli repulsion term +c use_xrepel logical flag governing use of exchange repulsion term c use_disp logical flag governing use of dispersion potential c use_charge logical flag governing use of charge-charge potential c use_chgdpl logical flag governing use of charge-dipole potential @@ -56,6 +57,7 @@ module potent logical use_pitors,use_strtor logical use_angtor,use_tortor logical use_vdw,use_repel + logical use_xrepel logical use_disp,use_charge logical use_chgdpl,use_dipole logical use_mpole,use_polar diff --git a/source/prmkey.f b/source/prmkey.f index dd81b90b..156efb64 100644 --- a/source/prmkey.f +++ b/source/prmkey.f @@ -132,7 +132,11 @@ subroutine prmkey (text) call getword (record,value,next) if (value .eq. 'ONLY') call potoff use_repel = .true. - if (value .eq. 'NONE') use_repel = .false. + use_xrepel = .true. + if (value .eq. 'NONE') then + use_repel = .false. + use_xrepel = .false. + end if else if (keyword(1:15) .eq. 'DISPERSIONTERM ') then call getword (record,value,next) if (value .eq. 'ONLY') call potoff @@ -549,6 +553,7 @@ subroutine potoff use_tortor = .false. use_vdw = .false. use_repel = .false. + use_xrepel = .false. use_disp = .false. use_charge = .false. use_chgdpl = .false. @@ -622,6 +627,7 @@ subroutine nbondoff c use_vdw = .false. use_repel = .false. + use_xrepel = .false. use_disp = .false. use_charge = .false. use_chgdpl = .false. diff --git a/source/readprm.f b/source/readprm.f index ac2cef91..cb9ef456 100644 --- a/source/readprm.f +++ b/source/readprm.f @@ -50,6 +50,7 @@ subroutine readprm use kurybr use kvdws use kvdwpr + use kxrepl use merck use params implicit none @@ -78,7 +79,7 @@ subroutine readprm integer tkey(maxtgrd2) real*8 wght,rd real*8 ep,rdn - real*8 spr,apr,epr + real*8 spr,apr,epr,cpr real*8 cdp,adp real*8 an1,an2,an3 real*8 ba1,ba2 @@ -2032,6 +2033,22 @@ subroutine readprm else if (ie .eq. 1) then mmffaroma(ia,if) = ic end if +c +c exchange repulsion parameters +c + else if (keyword(1:11) .eq. 'XREPULSION ') then + ia = 0 + spr = 0.0d0 + apr = 0.0d0 + cpr = 1.0d0 + string = record(next:240) + read (string,*,err=800,end=800) ia,spr,apr,cpr + 800 continue + if (ia .ne. 0) then + pxrz(ia) = spr + pxrdmp(ia) = apr + pxrcr(ia) = cpr + end if end if end do return diff --git a/source/respa.f b/source/respa.f index 6c71f483..5f8f8ba0 100644 --- a/source/respa.f +++ b/source/respa.f @@ -244,12 +244,14 @@ subroutine gradfast (energy,derivs) logical save_mpole,save_polar logical save_chgtrn,save_rxnfld logical save_solv,save_list + logical save_xrepel c c c save the original state of slow-evolving potentials c save_vdw = use_vdw save_repel = use_repel + save_xrepel = use_xrepel save_disp = use_disp save_charge = use_charge save_chgdpl = use_chgdpl @@ -265,6 +267,7 @@ subroutine gradfast (energy,derivs) c use_vdw = .false. use_repel = .false. + use_xrepel = .false. use_disp = .false. use_charge = .false. use_chgdpl = .false. @@ -284,6 +287,7 @@ subroutine gradfast (energy,derivs) c use_vdw = save_vdw use_repel = save_repel + use_xrepel = save_xrepel use_disp = save_disp use_charge = save_charge use_chgdpl = save_chgdpl diff --git a/source/testpair.f b/source/testpair.f index e0ac4a81..f7d59dce 100644 --- a/source/testpair.f +++ b/source/testpair.f @@ -89,7 +89,11 @@ program testpair npair = 0 nterm = 0 if (use_vdw) nterm = nterm + 1 - if (use_repel) nterm = nterm + 1 + if (use_repel) then + nterm = nterm + 1 + else if (use_xrepel) then + nterm = nterm + 1 + end if if (use_disp) nterm = nterm + 1 if (use_charge) nterm = nterm + 1 if (use_chgdpl) nterm = nterm + 1 @@ -269,7 +273,11 @@ program testpair if (vdwtyp .eq. 'BUFFERED-14-7') call ehal if (vdwtyp .eq. 'GAUSSIAN') call egauss end if - if (use_repel) call erepel + if (use_repel) then + call erepel + else if (use_xrepel) then + call exrepel + end if if (use_disp) call edisp if (use_charge) call echarge if (use_chgdpl) call echgdpl @@ -307,7 +315,11 @@ program testpair if (vdwtyp .eq. 'BUFFERED-14-7') call ehal if (vdwtyp .eq. 'GAUSSIAN') call egauss end if - if (use_repel) call erepel + if (use_repel) then + call erepel + else if (use_xrepel) then + call exrepel + end if if (use_disp) call edisp if (use_charge) call echarge if (use_chgdpl) call echgdpl @@ -342,7 +354,11 @@ program testpair if (vdwtyp .eq. 'BUFFERED-14-7') call ehal if (vdwtyp .eq. 'GAUSSIAN') call egauss end if - if (use_repel) call erepel + if (use_repel) then + call erepel + else if (use_xrepel) then + call exrepel + end if if (use_disp) call edisp if (use_charge) call echarge if (use_chgdpl) call echgdpl @@ -388,7 +404,11 @@ program testpair if (vdwtyp .eq. 'BUFFERED-14-7') call ehal1 if (vdwtyp .eq. 'GAUSSIAN') call egauss1 end if - if (use_repel) call erepel1 + if (use_repel) then + call erepel1 + else if (use_xrepel) then + call exrepel1 + end if if (use_disp) call edisp1 if (use_charge) call echarge1 if (use_chgdpl) call echgdpl1 @@ -442,7 +462,11 @@ program testpair if (vdwtyp .eq. 'BUFFERED-14-7') call ehal1 if (vdwtyp .eq. 'GAUSSIAN') call egauss1 end if - if (use_repel) call erepel1 + if (use_repel) then + call erepel1 + else if (use_xrepel) then + call exrepel1 + end if if (use_disp) call edisp1 if (use_charge) call echarge1 if (use_chgdpl) call echgdpl1 @@ -493,7 +517,11 @@ program testpair if (vdwtyp .eq. 'BUFFERED-14-7') call ehal1 if (vdwtyp .eq. 'GAUSSIAN') call egauss1 end if - if (use_repel) call erepel1 + if (use_repel) then + call erepel1 + else if (use_xrepel) then + call exrepel1 + end if if (use_disp) call edisp1 if (use_charge) call echarge1 if (use_chgdpl) call echgdpl1 diff --git a/source/xrepel.f b/source/xrepel.f new file mode 100644 index 00000000..7d747abf --- /dev/null +++ b/source/xrepel.f @@ -0,0 +1,38 @@ +c +c +c ############################################################ +c ## COPYRIGHT (C) 2022 by Moses KJ Chung & Jay W. Ponder ## +c ## All Rights Reserved ## +c ############################################################ +c +c ############################################################### +c ## ## +c ## module xrepel -- exch repulsion for current structure ## +c ## ## +c ############################################################### +c +c +c nxrep total number of repulsion sites in the system +c ixrep number of the atom for each repulsion site +c xreplist repulsion multipole site for each atom (0=none) +c zpxr nuclear charge parameter value for each atom +c dmppxr exchange repulsion alpha damping value for each atom +c crpxr ratio of p/s orbital cofficients for each atom +c cpxr local pseudo wavefunction coefficient for each atom +c rcpxr global pseudo wavefunction coefficient for each atom +c xrepole repulsion Cartesian multipoles in the local frame +c +c + module xrepel + implicit none + integer nxrep + integer, allocatable :: ixrep(:) + integer, allocatable :: xreplist(:) + real*8, allocatable :: zpxr(:) + real*8, allocatable :: dmppxr(:) + real*8, allocatable :: crpxr(:) + real*8, allocatable :: cpxr(:,:) + real*8, allocatable :: rcpxr(:,:) + real*8, allocatable :: xrepole(:,:) + save + end