diff --git src/pmemd.amba/src/amoeba_direct.fpp src/pmemd.amba/src/amoeba_direct.fpp index 7ab36dd..cff3497 100644 --- src/pmemd.amba/src/amoeba_direct.fpp +++ src/pmemd.amba/src/amoeba_direct.fpp @@ -2073,7 +2073,8 @@ subroutine am_vdw_direct_ene_frc_i(atm_i, img, ipairs_sublst, x_tran, & pair_cnt, crd, ene_vdw, frc, img_frc, & virial, img_atm_map) - use mdin_amoeba_dat_mod, only : do_vdw_taper + use mdin_amoeba_dat_mod, only : do_vdw_taper, soft_lambda, soft_expo, soft_alpha, dipole_scf_iter_max, & + soft_atom_range1,soft_atom_range2, soft_line use amoeba_flags_mod use img_mod @@ -2098,6 +2099,7 @@ subroutine am_vdw_direct_ene_frc_i(atm_i, img, ipairs_sublst, x_tran, & integer :: itran, it, jt, ih, jh, idx integer :: sublst_idx integer :: atm_j, img_j + integer :: i_range double precision :: wi, wj double precision :: delx, dely, delz double precision :: delr, delr2 @@ -2193,6 +2195,16 @@ subroutine am_vdw_direct_ene_frc_i(atm_i, img, ipairs_sublst, x_tran, & dt1drho = -7.d0 * t1 / (rho + vdw_buf_delta) dt2drho = -7.d0 * t2 * (rho6 / (rho7 + vdw_buf_gamma)) drhodr = 1.d0 / rad + do i_range = 1,soft_line + if(((atm_i.ge.soft_atom_range1(i_range)).and.(atm_i.le.soft_atom_range2(i_range))) & + .or. ((atm_j.ge.soft_atom_range1(i_range)) .and.(atm_j.le.soft_atom_range2(i_range)))) then + eps = eps * soft_lambda ** soft_expo + t1 = (1.d0 + vdw_buf_delta)**7 / (soft_alpha * (1.d0 - soft_lambda)**2 + (rho + vdw_buf_delta)**7) + t2 = (1.d0 + vdw_buf_gamma) / (soft_alpha * (1.d0 - soft_lambda)**2 + rho**7 + vdw_buf_gamma) + dt1drho = (-7.d0 * (rho + vdw_buf_delta)**6 * t1) / (soft_alpha *(1.d0 - soft_lambda)**2 + (rho + vdw_buf_delta)**7) + dt2drho = (-7.d0 * rho**6 * t2) / (soft_alpha * (1.d0 - soft_lambda)**2 + rho**7 + vdw_buf_gamma) + endif + end do f = eps * t1 * (t2 - 2.d0) dfdr = eps * (dt1drho * (t2 - 2.d0) + t1 * dt2drho) * drhodr diff --git src/pmemd.amba/src/amoeba_induced.fpp src/pmemd.amba/src/amoeba_induced.fpp index 1bde607..5896025 100644 --- src/pmemd.amba/src/amoeba_induced.fpp +++ src/pmemd.amba/src/amoeba_induced.fpp @@ -285,6 +285,18 @@ subroutine amoeba_induced_dat_final_setup(num_ints, num_reals) if (polarizability(n) .gt. polmin) then is_polarizable(n) = .true. atm_qterm(n) = 1.d0 / sqrt(polarizability(n)) + if ( abs(atm_mass(n) - 40.078) < 0.0001d0 ) & + write (6,*) "Ca damp",n,atm_mass(n),polarizability(n) + if ( abs(atm_mass(n) - 40.078) < 0.0001d0 ) & + atm_qterm(n) = 1 /1.35d0**3*atm_qterm(n) + if ( abs(atm_mass(n) - 24.305) < 0.0001d0 ) & + write (6,*) "Mg damp",n,atm_mass(n),polarizability(n) + if ( abs(atm_mass(n) - 24.305) < 0.0001d0 ) & + atm_qterm(n) = 1 /1.60d0**3*atm_qterm(n) + if ( abs(atm_mass(n) - 65.409) < 0.0001d0 ) & + write (6,*) "Zn damp",n,atm_mass(n),polarizability(n) + if ( abs(atm_mass(n) - 65.409) < 0.0001d0 ) & + atm_qterm(n) = 1 /1.23d0**3*atm_qterm(n) else is_polarizable(n) = .false. atm_qterm(n) = 1.d0 / sqrt(polmin) diff --git src/pmemd.amba/src/amoeba_vdw.fpp src/pmemd.amba/src/amoeba_vdw.fpp index 3710051..99d0f16 100644 --- src/pmemd.amba/src/amoeba_vdw.fpp +++ src/pmemd.amba/src/amoeba_vdw.fpp @@ -369,7 +369,9 @@ end subroutine am_vdw_switch subroutine am_vdw_longrange_factor(atm_cnt) - use mdin_amoeba_dat_mod, only : do_vdw_taper + use mdin_amoeba_dat_mod, only : do_vdw_taper, soft_lambda,AMOEBA_read_soft, & + soft_atom_range1,soft_atom_range2,soft_line,& + vdw_longrange_lambda use gbl_constants_mod, only : PI implicit none @@ -381,17 +383,29 @@ subroutine am_vdw_longrange_factor(atm_cnt) ! Local variables: integer :: ier, n, nt, kdel, ndel, i, j + integer :: i_range, i_soft, i_soft_type double precision :: r, r1, r2, r3, r4, r5 double precision :: req, eps, f, sume, sumv, t1, t2, rho, switch, delr double precision :: dt1drho, dt2drho, drhodr, dfdr, dswitch_dr, g1, g2 + integer,save,allocatable :: lig_type_ct(:) + if (vdw_longrange_lambda .ne. 1.0 .or. soft_lambda .ne. 1.0) then + call AMOEBA_read_soft() + endif + allocate(lig_type_ct(vdw_param_cnt),stat=ier) do n = 1, vdw_param_cnt vdw_type_count(n) = 0 + lig_type_ct(n) = 0 end do do n = 1, atm_cnt nt = vdw_atom_type(n) vdw_type_count(nt) = vdw_type_count(nt) + 1 + do i_range = 1, soft_line + if (n .ge. soft_atom_range1(i_range).and. n.le.soft_atom_range2(i_range)) then + lig_type_ct(nt) = lig_type_ct(nt) + 1 + endif + end do end do ! Choose r2 to be 60 Angstroms or ~20 sigma. @@ -458,11 +472,12 @@ subroutine am_vdw_longrange_factor(atm_cnt) ! Note the 2 * pi below not 4 * pi---since we do each i, j pair 2x. - ene_vdw_longrange_factor = ene_vdw_longrange_factor + 2.d0 * PI * & - vdw_type_count(i) * vdw_type_count(j) * sume + ene_vdw_longrange_factor = ene_vdw_longrange_factor + 2.d0 * PI * sume * (vdw_type_count(i)*vdw_type_count(j) & + -(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j))) + + vir_vdw_longrange_factor = vir_vdw_longrange_factor + 2.d0 * PI * sumv * (vdw_type_count(i)*vdw_type_count(j) & + -(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j))) - vir_vdw_longrange_factor = vir_vdw_longrange_factor + 2.d0 * PI * & - vdw_type_count(i) * vdw_type_count(j) * sumv end do end do diff --git src/pmemd.amba/src/mdin_amoeba_dat.fpp src/pmemd.amba/src/mdin_amoeba_dat.fpp index f2bc655..cc43126 100644 --- src/pmemd.amba/src/mdin_amoeba_dat.fpp +++ src/pmemd.amba/src/mdin_amoeba_dat.fpp @@ -14,7 +14,7 @@ module mdin_amoeba_dat_mod ! Data that should be broadcast to slave processes from the master: - integer, parameter :: mdin_amoeba_int_cnt = 23 + integer, parameter :: mdin_amoeba_int_cnt = 24 integer do_amoeba_valence, do_amoeba_nonbond, do_bond, & do_ureyb, do_reg_angle, do_trig_angle, do_opbend, & @@ -23,7 +23,9 @@ module mdin_amoeba_dat_mod do_adjust, do_direct, do_self, do_vdw, & do_induced, do_vdw_taper, do_vdw_longrange, & beeman_integrator, dipole_scf_iter_max, & - amoeba_verbose + amoeba_verbose, soft_expo + + integer :: soft_atom_range1(20), soft_atom_range2(20), soft_line common / mdin_amoeba_int / do_amoeba_valence, do_amoeba_nonbond, do_bond, & do_ureyb, do_reg_angle, do_trig_angle, do_opbend, & @@ -32,20 +34,23 @@ module mdin_amoeba_dat_mod do_adjust, do_direct, do_self, do_vdw, & do_induced, do_vdw_taper, do_vdw_longrange, & beeman_integrator, dipole_scf_iter_max, & - amoeba_verbose + amoeba_verbose, soft_expo + common soft_atom_range1, soft_atom_range2, soft_line save :: / mdin_amoeba_int / - integer, parameter :: mdin_amoeba_dbl_cnt = 7 + integer, parameter :: mdin_amoeba_dbl_cnt = 10 double precision compress, dipole_scf_tol, ee_dsum_cut, & ee_damped_cut, sor_coefficient, & - thole_expon_coeff, vdw_taper + thole_expon_coeff, vdw_taper, & + soft_lambda, soft_alpha, vdw_longrange_lambda common / mdin_amoeba_dbl / compress, dipole_scf_tol, ee_dsum_cut, & ee_damped_cut, sor_coefficient, & - thole_expon_coeff, vdw_taper + thole_expon_coeff, vdw_taper, & + soft_lambda, soft_alpha, vdw_longrange_lambda save :: / mdin_amoeba_dbl / @@ -74,11 +79,12 @@ subroutine init_mdin_amoeba_dat do_trig_angle, do_opbend, do_torsion, do_str_torsion, & do_pi_torsion, do_strbend, do_torsion_torsion, & do_recip, do_adjust, do_direct, do_self, & - do_vdw, do_induced, amoeba_verbose, beeman_integrator, & + do_vdw, do_induced, amoeba_verbose, beeman_integrator, & sor_coefficient, dipole_scf_tol, & dipole_scf_iter_max, ee_dsum_cut, ee_damped_cut, & thole_expon_coeff, vdw_taper, do_vdw_taper, & - do_vdw_longrange, compress + do_vdw_longrange, compress, & + soft_lambda, soft_alpha, soft_expo, vdw_longrange_lambda ! Set default values. @@ -113,6 +119,10 @@ subroutine init_mdin_amoeba_dat sor_coefficient = 0.75d0 thole_expon_coeff = 0.39d0 vdw_taper = 0.9d0 + soft_lambda = 1.d0 + soft_alpha = 0.5d0 + soft_expo = 4 + vdw_longrange_lambda = 1.d0 rewind(mdin) call nmlsrc('amoeba', mdin, ifind) @@ -156,6 +166,47 @@ subroutine bcast_mdin_amoeba_dat end subroutine bcast_mdin_amoeba_dat #endif +subroutine AMOEBA_read_soft + + use parallel_dat_mod + + integer pos_dash,length,i_range,i + character(len=20) atm_range + character(len=6) temp + logical alive + inquire(file='soft_atm.txt',exist=alive) + if (alive) then + do i=1,20 + soft_atom_range1(i)=0; + soft_atom_range2(i)=0; + end do + soft_line=0 + + open (11,FILE='soft_atm.txt') + do while (.true.) + soft_line=soft_line+1 + read(11,*,end=99) atm_range + pos_dash=scan(atm_range,'-') + length=len_trim(atm_range) + if(pos_dash.gt.0) then + temp=atm_range(1:pos_dash-1) + read(temp,*) soft_atom_range1(soft_line) + temp=atm_range(pos_dash+1:length) + read(temp,*) soft_atom_range2(soft_line) + else + read(atm_range,*) soft_atom_range1(soft_line) + read(atm_range,*) soft_atom_range2(soft_line) + endif + end do + else + write(6,*)'soft_atm.txt not found' + call mexit(6, 1) + endif + 99 continue + close(11) + +end subroutine AMOEBA_read_soft + !******************************************************************************* ! ! Subroutine: validate_mdin_amoeba_dat diff --git src/sander/amoeba_induced.f src/sander/amoeba_induced.f index eb3e3d8..d0d1aba 100644 --- src/sander/amoeba_induced.f +++ src/sander/amoeba_induced.f @@ -110,13 +110,14 @@ subroutine AM_INDUCED_bcast(num_atoms) end subroutine AM_INDUCED_bcast #endif -function AM_INDUCED_readparm(nf,num_atoms) +function AM_INDUCED_readparm(nf,num_atoms,mass) integer :: AM_INDUCED_readparm integer,intent(in) :: nf,num_atoms #include "do_flag.h" integer :: j,n,ier,dim1,numpolar _REAL_ :: polmin + _REAL_ :: mass(*) AM_INDUCED_readparm = 0 polmin = 1.d-12 @@ -177,7 +178,19 @@ function AM_INDUCED_readparm(nf,num_atoms) if ( polarizability(n) > polmin )then is_polarizable(n) = .true. sq_polinv(n) = 1.d0 / sqrt(polarizability(n)) - endif + if ( abs(mass(n) - 40.078) < 0.0001d0 ) & + write (6,*) "Ca damp",n,mass(n),polarizability(n) + if ( abs(mass(n) - 40.078) < 0.0001d0 ) & + sq_polinv(n) = 1 /1.35d0**3*sq_polinv(n) + if ( abs(mass(n) - 24.305) < 0.0001d0 ) & + write (6,*) "Mg damp",n,mass(n),polarizability(n) + if ( abs(mass(n) - 24.305) < 0.0001d0 ) & + sq_polinv(n) = 1 /1.60d0**3*sq_polinv(n) + if ( abs(mass(n) - 65.409) < 0.0001d0 ) & + write (6,*) "Zn damp",n,mass(n),polarizability(n) + if ( abs(mass(n) - 65.409) < 0.0001d0 ) & + sq_polinv(n) = 1 /1.34d0**3*sq_polinv(n) + endif enddo AM_INDUCED_readparm = 1 do_flag = ibset(do_flag,VALID_BIT) diff --git src/sander/amoeba_interface.f src/sander/amoeba_interface.f index 283c502..286f57e 100644 --- src/sander/amoeba_interface.f +++ src/sander/amoeba_interface.f @@ -8,7 +8,7 @@ module amoeba_interface AMOEBA_deallocate,AM_VAL_eval,AM_NonBond_eval contains !------------------------------------------------------------------------------- -subroutine AM_NONBOND_readparm(nf,natom) +subroutine AM_NONBOND_readparm(nf,natom,mass) use amoeba_multipoles, only : AM_MPOLE_readparm use amoeba_adjust, only : AM_ADJUST_readparm use amoeba_vdw, only : AM_VDW_read_parm @@ -18,10 +18,11 @@ subroutine AM_NONBOND_readparm(nf,natom) integer,intent(in) :: nf,natom integer :: mpole_valid,adjust_valid,vdw_valid,polar_valid + _REAL_ mass(*) mpole_valid = AM_MPOLE_readparm(nf,natom) adjust_valid = AM_ADJUST_readparm(nf) vdw_valid = AM_VDW_read_parm(nf) - polar_valid = AM_INDUCED_readparm(nf,natom) + polar_valid = AM_INDUCED_readparm(nf,natom,mass) ! alocate to initialize recip call AM_RECIP_allocate(natom) call AM_RUNMD_init(natom) @@ -65,15 +66,16 @@ subroutine AM_VAL_readparm(nf,ntf,ntc) endif end subroutine AM_VAL_readparm !------------------------------------------------------------------------------- -subroutine AMOEBA_readparm(nf,ntf,ntc,numatoms) +subroutine AMOEBA_readparm(nf,ntf,ntc,numatoms,mass) use amoeba_mdin, only : iamoeba integer,intent(in) :: nf,numatoms integer,intent(inout) :: ntf,ntc !set these to 8,1 if amoeba-friendly prmtop !so regular amber valence & shake not run + _REAL_ mass(*) call AMOEBA_check_parm_legal(nf) if ( iamoeba /= 1 )return call AM_VAL_readparm(nf,ntf,ntc) - call AM_NONBOND_readparm(nf,numatoms) + call AM_NONBOND_readparm(nf,numatoms,mass) end subroutine AMOEBA_readparm !------------------------------------------------------------------------------- subroutine AM_VAL_deallocate() diff --git src/sander/amoeba_mdin.f src/sander/amoeba_mdin.f index 9d9f217..f8e4247 100644 --- src/sander/amoeba_mdin.f +++ src/sander/amoeba_mdin.f @@ -18,8 +18,16 @@ module amoeba_mdin _REAL_,save :: vdw_taper = 0.9d0 _REAL_,save :: thole_expon_coeff=0.39d0 _REAL_,save :: compress = 0.000046d0 + _REAL_,save :: soft_lambda = 1.0d0 + _REAL_,save :: soft_alpha = 0.5d0 + _REAL_,save :: soft_expo = 4 + _REAL_,save :: vdw_longrange_lambda = 1.0d0 integer, save :: amoeba_verbose = 0 integer,save :: beeman_integrator = 0 + + !variables and arrays for softcore file + integer,save :: soft_atom_range1(20),soft_atom_range2(20),soft_line + public AMOEBA_read_mdin,iamoeba,do_amoeba_valence, & do_amoeba_nonbond,amoeba_verbose,beeman_integrator, & sor_coefficient,dipole_scf_tol,dipole_scf_iter_max, & @@ -27,7 +35,10 @@ module amoeba_mdin dipole_scf_use_cg, dipole_scf_cg_niter, & #endif /* DISABLE_AMOEBA_CG */ ee_dsum_cut,ee_damped_cut,thole_expon_coeff,vdw_taper, & - compress,do_vdw_taper,do_vdw_longrange,am_nbead + compress,do_vdw_taper,do_vdw_longrange,am_nbead, & + soft_lambda,soft_alpha,soft_expo, & + AMOEBA_read_soft, soft_atom_range1,soft_atom_range2, soft_line, & + vdw_longrange_lambda contains !------------------------------------------------------------------------------- subroutine AMOEBA_read_mdin(nf) @@ -53,7 +64,8 @@ subroutine AMOEBA_read_mdin(nf) #endif /* DISABLE_AMOEBA_CG */ ee_dsum_cut,ee_damped_cut, & thole_expon_coeff,vdw_taper,do_vdw_taper,do_vdw_longrange, & - compress,am_nbead + compress,am_nbead,& + soft_lambda,soft_alpha,soft_expo,vdw_longrange_lambda # include "files.h" @@ -67,5 +79,52 @@ subroutine AMOEBA_read_mdin(nf) end subroutine AMOEBA_read_mdin !---------------------------------------------------------- +subroutine AMOEBA_read_soft() + integer pos_dash,length,i_range,i + character(len=20) atm_range + character(len=6) temp + + !flag whether soft_atm.txt exists + logical alive + inquire(file='soft_atm.txt',exist=alive) + if(alive) then + do i=1,20 + soft_atom_range1(i)=0; + soft_atom_range2(i)=0; + end do + + soft_line=0 + open (11,FILE='soft_atm.txt') + do while (.true.) + + read(11,*,end=99) atm_range + !write(*,*) atm_range + pos_dash=scan(atm_range,'-') + !write(*,*) "position",pos_dash + length=len_trim(atm_range) + if (length > 0) then + soft_line=soft_line+1 + endif + !write(*,*) "length",length + if(pos_dash.gt.0) then + temp=atm_range(1:pos_dash-1) + read(temp,*) soft_atom_range1(soft_line) + temp=atm_range(pos_dash+1:length) + read(temp,*) soft_atom_range2(soft_line) + + else + read(atm_range,*) soft_atom_range1(soft_line) + read(atm_range,*) soft_atom_range2(soft_line) + endif + end do + else + write(6,*) 'soft_atm.txt not found' + call mexit(6,1) + endif + 99 continue + close(11) + +end subroutine AMOEBA_read_soft +!--------------------------------------------------------- end module amoeba_mdin !------------------------------------------------------------------------------- diff --git src/sander/amoeba_vdw.f src/sander/amoeba_vdw.f index 991d6c7..575c907 100644 --- src/sander/amoeba_vdw.f +++ src/sander/amoeba_vdw.f @@ -178,21 +178,36 @@ use amoeba_mdin, only : vdw_taper end subroutine AM_VDW_switch !------------------------------------------------------------ subroutine AM_VDW_longrange_factor(num_atoms) - use amoeba_mdin, only : do_vdw_taper + use amoeba_mdin, only : do_vdw_taper, & + soft_lambda,soft_atom_range1,soft_atom_range2,& + soft_line,vdw_longrange_lambda,AMOEBA_read_soft use constants, only : zero,one,two,three,four,five,seven,pi integer,intent(in) :: num_atoms - integer ier,n,nt,kdel,ndel,i,j + integer ier,n,nt,kdel,ndel,i,j,i_range,i_soft,i_soft_type + integer,save,allocatable :: lig_type_ct(:) _REAL_ :: r,r1,r2,r3,r4,r5,req,eps,f,sume,sumv,t1,t2,rho,switch,delr _REAL_ :: dt1drho,dt2drho,drhodr,dfdr,dswitch_dr,g1,g2 allocate(vdw_type_count(num_vdw_atom_types),stat=ier) REQUIRE(ier==0) + allocate(lig_type_ct(num_vdw_atom_types),stat=ier) + REQUIRE(ier==0) + + if (vdw_longrange_lambda .ne. 1.0 .or. soft_lambda .ne. 1.0) then + call AMOEBA_read_soft() + endif do n = 1,num_vdw_atom_types vdw_type_count(n) = 0 + lig_type_ct(n) = 0 enddo do n = 1,num_atoms nt = vdw_atom_type(n) vdw_type_count(nt) = vdw_type_count(nt) + 1 + do i_range = 1, soft_line + if (n .ge. soft_atom_range1(i_range).and. n.le.soft_atom_range2(i_range)) then + lig_type_ct(nt) = lig_type_ct(nt) + 1 + endif + enddo enddo ! choose r2 to be 60 Angstroms or ~20 sigma ! choose delr to be 0.05 angstroms, or ndel = 20*(r2-r1) @@ -247,10 +262,10 @@ subroutine AM_VDW_longrange_factor(num_atoms) endif enddo ! note the 2*pi below not 4*pi---since we do each i,j pair 2x - ene_vdw_longrange_factor = ene_vdw_longrange_factor + two*pi* & - vdw_type_count(i)*vdw_type_count(j)*sume - vir_vdw_longrange_factor = vir_vdw_longrange_factor + two*pi* & - vdw_type_count(i)*vdw_type_count(j)*sumv + ene_vdw_longrange_factor = ene_vdw_longrange_factor + two*pi*sume*(vdw_type_count(i)*vdw_type_count(j) & + -(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j))) + vir_vdw_longrange_factor = vir_vdw_longrange_factor + two*pi*sumv*(vdw_type_count(i)*vdw_type_count(j) & + -(1-vdw_longrange_lambda)*(lig_type_ct(i)*vdw_type_count(j)+(vdw_type_count(i)-lig_type_ct(i))*lig_type_ct(j))) enddo enddo @@ -384,12 +399,14 @@ subroutine AM_VDW_DIRECT_ene_frc_i(i,ipairs,numtot,xk,yk,zk, & crd,ene_vdw,frc,virial) use nblist, only: bckptr,imagcrds,tranvec use constants, only : zero,one,two,three,four,five,seven - use amoeba_mdin, only : do_vdw_taper + use amoeba_mdin, only : do_vdw_taper, & + soft_lambda,soft_alpha,soft_expo, & + soft_atom_range1,soft_atom_range2, soft_line integer,intent(in) :: i,ipairs(*),numtot _REAL_,intent(in) :: xk,yk,zk,crd(3,*) _REAL_,intent(inout) :: ene_vdw,frc(3,*),virial(3,3) - integer :: itran,it,jt,ih,jh,j,m,np,mask27 + integer :: itran,it,jt,ih,jh,j,m,np,mask27,i_range _REAL_ :: xktran(3,18) _REAL_ :: wi,wj,xi,yi,zi,delx,dely,delz,delr2,delr,eps,rad _REAL_ :: rho,t1,t2,dt1drho,dt2drho,drhodr,term,dfx,dfy,dfz @@ -448,6 +465,22 @@ subroutine AM_VDW_DIRECT_ene_frc_i(i,ipairs,numtot,xk,yk,zk, & dt1drho = -seven*t1 / (rho + vdw_buf_delta) dt2drho = -seven*t2 * (rho6 / (rho7 + vdw_buf_gamma)) drhodr = one / rad + do i_range = 1,soft_line +! if ((i.ge.soft_atom_range1(i_range).and. & +! i.le.soft_atom_range2(i_range) ) & +! .xor. & +! (j.ge.soft_atom_range1(i_range).and. & +! j.le.soft_atom_range2(i_range)) )then + if ( xor((i.ge.soft_atom_range1(i_range).and.i.le.soft_atom_range2(i_range)), & + (j.ge.soft_atom_range1(i_range).and.j.le.soft_atom_range2(i_range)) ) )then + eps = eps * soft_lambda ** soft_expo + t1 = (one + vdw_buf_delta)**7 / (soft_alpha * (1 - soft_lambda)**2 + (rho + vdw_buf_delta)**7) + t2 = (one + vdw_buf_gamma) / (soft_alpha * (1 - soft_lambda)**2 + rho7 + vdw_buf_gamma) + dt1drho = (-seven * (rho + vdw_buf_delta)**6 * t1) / (soft_alpha * (1 - soft_lambda)**2 & + + (rho + vdw_buf_delta)**7) + dt2drho = (-seven * rho6 * t2) / (soft_alpha * (1 - soft_lambda)**2 + rho7 + vdw_buf_gamma) + endif + end do f = eps*t1*(t2 - two) dfdr = eps*(dt1drho*(t2 - two) + t1*dt2drho)*drhodr if ( do_vdw_taper == 1 .and. delr2 > vdw_switch_on_2 )then diff --git src/sander/sander.f src/sander/sander.f index fafe8a8..bd92276 100644 --- src/sander/sander.f +++ src/sander/sander.f @@ -316,7 +316,7 @@ subroutine sander() ! --- finish reading the prmtop file and other user input: call rdparm2(x,ix,ih,ipairs,8) - call AMOEBA_readparm(8,ntf,ntc,natom)! ntf,ntc get reset if amoeba prmtop + call AMOEBA_readparm(8,ntf,ntc,natom,x(lmass))! ntf,ntc get reset if amoeba prmtop #ifdef _XRAY call xray_read_parm(8,6) #endif