From 71b2f738d477dfd7e459fafdb5444a5038dfb1cd Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Fri, 2 Jun 2006 17:27:06 +0000 Subject: [PATCH] Bug fix for onfo. Fully conforming amber results to FIST results. introduced few other tests to check that amber and fist will be compatible in the future.. Introduced a keyword to avoid the scaling of the cutoff energy.. In this way amber and fist give exactly the same results.A bug fix also in the topology_converter.. cosmetics around.. reset 2 regtests. svn-origin-rev: 5204 --- src/fist_intra_force.F | 9 +- src/force_field_types.F | 1 + src/force_fields.F | 188 +++++++++++----------- src/input_cp2k_mm.F | 7 + src/pair_potential.F | 36 +++-- src/splines.F | 15 +- src/topology_util.F | 40 +---- tests/Fist/regtest/TEST_FILES | 2 + tests/Fist/regtest/TEST_FILES_RESET | 4 + tests/Fist/regtest/gly.psf | 83 ++++++++++ tests/Fist/regtest/gly_amber.inp | 61 +++++++ tests/Fist/regtest/h2o2_amber.inp | 5 +- tests/Fist/regtest/mol.psf | 36 +++++ tests/Fist/regtest/mol_amber.inp | 80 +++++++++ tests/Fist/sample_pdb/gly.pdb | 22 +++ tests/Fist/sample_pdb/mol.pdb | 9 ++ tests/Fist/sample_pot/gly.pot | 165 +++++++++++++++++++ tests/Fist/sample_pot/mol.pot | 94 +++++++++++ tools/Fist/topology_converter/leap2fist.F | 112 ++++++++++++- 19 files changed, 815 insertions(+), 154 deletions(-) create mode 100644 tests/Fist/regtest/gly.psf create mode 100644 tests/Fist/regtest/gly_amber.inp create mode 100644 tests/Fist/regtest/mol.psf create mode 100644 tests/Fist/regtest/mol_amber.inp create mode 100644 tests/Fist/sample_pdb/gly.pdb create mode 100644 tests/Fist/sample_pdb/mol.pdb create mode 100644 tests/Fist/sample_pot/gly.pot create mode 100644 tests/Fist/sample_pot/mol.pot diff --git a/src/fist_intra_force.F b/src/fist_intra_force.F index 26b5528dcf..b31ecda169 100644 --- a/src/fist_intra_force.F +++ b/src/fist_intra_force.F @@ -110,6 +110,7 @@ SUBROUTINE force_intra_control ( molecule_set, molecule_kind_set, local_molecule nimptors, nkind, nmol_per_kind, nonefs, ntorsions, nub REAL(KIND=dp) :: d12, d32, dist, dist1, dist2, energy, fscalar, id12, & id32, is32, ism, isn, iw1, iw2, rab2, s32, sm, sn, theta, w1, w2 + REAL(KIND=dp) :: rinv, rinv2 REAL(KIND=dp), DIMENSION(3) :: b12, b32, g1, g2, g3, gt1, & gt2, gt3, gt4, rij, t12, t32, & t34, t43, tm, tn, u1, u2, & @@ -279,7 +280,7 @@ SUBROUTINE force_intra_control ( molecule_set, molecule_kind_set, local_molecule index_a = torsion_list ( itorsion ) % a + first_atom - 1 index_b = torsion_list ( itorsion ) % b + first_atom - 1 index_c = torsion_list ( itorsion ) % c + first_atom - 1 - index_d = torsion_list ( itorsion ) % d + first_atom - 1 + index_d = torsion_list ( itorsion ) % d + first_atom - 1 t12 = particle_set ( index_a ) % r - particle_set ( index_b ) % r t32 = particle_set ( index_c ) % r - particle_set ( index_b ) % r t43 = particle_set ( index_d ) % r - particle_set ( index_c ) % r @@ -415,8 +416,14 @@ SUBROUTINE force_intra_control ( molecule_set, molecule_kind_set, local_molecule rij = particle_set ( index_b ) % r - particle_set ( index_a ) % r rab2 = DOT_PRODUCT ( rij, rij ) spline_data => onef_list ( ionef ) % onfo_kind % spline_data + ! LJ contribution energy = potential_s ( spline_data, rab2, fscalar ) pot_onef = pot_onef + energy + ! EI contribution + rinv2 = 1.0_dp/rab2 + rinv = SQRT(rinv2) + pot_onef = pot_onef + spline_data%ei_scale14 * rinv - spline_data%ei_cutoff14 + fscalar = fscalar + spline_data%ei_scale14 * rinv2 * rinv particle_set ( index_a )% f(1) = particle_set ( index_a ) % f(1) -rij(1) * fscalar particle_set ( index_a )% f(2) = particle_set ( index_a ) % f(2) -rij(2) * fscalar particle_set ( index_a )% f(3) = particle_set ( index_a ) % f(3) -rij(3) * fscalar diff --git a/src/force_field_types.F b/src/force_field_types.F index 4a2f983e56..1621fd2ffc 100644 --- a/src/force_field_types.F +++ b/src/force_field_types.F @@ -168,6 +168,7 @@ MODULE force_field_types max_energy,& emax_spline,& rlow_nb + LOGICAL :: scale_cutoff END TYPE force_field_type ! *** Public subroutines *** diff --git a/src/force_fields.F b/src/force_fields.F index ed9181bc0e..21369a6230 100644 --- a/src/force_fields.F +++ b/src/force_fields.F @@ -334,6 +334,7 @@ CONTAINS CALL section_vals_val_get(ff_section,"SPLINE%EPS_SPLINE",r_val=ff_type%eps_spline,error=error) CALL section_vals_val_get(ff_section,"SPLINE%EMAX_SPLINE",r_val=ff_type%emax_spline,error=error) CALL section_vals_val_get(ff_section,"SPLINE%EMAX_ACCURACY",r_val=ff_type%max_energy,error=error) + CALL section_vals_val_get(ff_section,"SCALE_CUTOFF",l_val=ff_type%scale_cutoff,error=error) tmp_section => section_vals_get_subs_vals(ff_section,"NONBONDED",error=error) CALL section_vals_get(tmp_section,explicit=explicit,error=error) @@ -2712,7 +2713,96 @@ CONTAINS CALL timestop(handle2) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 12. Set spline_bond_env, use for ONFO interaction + !12. Set up atomic_kind_set()%elp_potentail%qeff (PART 2) + !----------------------------------------------------------------------------- + CALL timeset(routineN//'_qeff',handle2) + ALLOCATE(charge(SIZE(atomic_kind_set)),STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + DO i=1,SIZE(atomic_kind_set) + atomic_kind => atomic_kind_set(i) + CALL get_atomic_kind(atomic_kind=atomic_kind,& + elp_potential=elp_potential,& + atom_list=my_atom_list,& + name=atmname) + CALL get_potential(potential=elp_potential,& + qeff=charge(i)) + + found = .FALSE. + only_qm = qmmm_ff_precond_only_qm(id1=atmname,is_link=is_link_atom) + CALL uppercase(atmname) + IF(charge(i)/=-HUGE(0.0_dp)) found = .TRUE. + ! loop over params from charmm + IF(ASSOCIATED(ff_type%chm_info%charge_atm)) THEN + DO j=1,SIZE(ff_type%chm_info%charge_atm) + IF(TRIM(ff_type%chm_info%charge_atm(j))==atmname) THEN + charge(i) = ff_type%chm_info%charge(j) + IF (found) CALL Pwarning("Multiple charge declarations: "//TRIM(atmname)//& + " overwriting!",globenv,error) + found = .TRUE. + END IF + END DO + END IF + ! always have the input param last to overwrite all the other ones + IF(ASSOCIATED(ff_type%inp_info%charge_atm)) THEN + DO j=1,SIZE(ff_type%inp_info%charge_atm) + IF (iw>0) WRITE(iw,*)"Charge Checking ::",TRIM(ff_type%inp_info%charge_atm(j)),atmname + IF(TRIM(ff_type%inp_info%charge_atm(j))==atmname) THEN + charge(i) = ff_type%inp_info%charge(j) + IF (found) CALL Pwarning("Multiple charge declarations: "//TRIM(atmname)//& + " overwriting!",globenv,error) + found = .TRUE. + END IF + END DO + END IF + IF (ASSOCIATED(ff_type%inp_info%nonbonded)) THEN + IF (ASSOCIATED(ff_type%inp_info%nonbonded%pot)) THEN + ! Let's find the nonbonded potential where this atom is involved + DO j = 1, SIZE(ff_type%inp_info%nonbonded%pot) + IF (atmname==ff_type%inp_info%nonbonded%pot(j)%pot %at1.OR.& + atmname==ff_type%inp_info%nonbonded%pot(j)%pot %at2) THEN + SELECT CASE(ff_type%inp_info%nonbonded%pot(j)%pot%type) + CASE (ea_type) + ! Nothing to complain.. we don't need charges + CASE default + IF(.NOT.found) THEN + CALL stop_program("FF_coordinate_pack","missing qeff parm "//atmname) + END IF + END SELECT + END IF + END DO + END IF + END IF + ! + ! QM/MM modifications + ! + IF (only_qm.AND.my_qmmm) THEN + IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN + scale_factor = 0.0_dp + IF (is_link_atom) THEN + ! + ! Find the scaling factor... + ! + DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) + IF (ANY(my_atom_list == qmmm_env%mm_link_atoms(ilink))) EXIT + END DO + CPPostcondition(ilink <= SIZE(qmmm_env%mm_link_atoms),cp_failure_level,routineP,error,failure) + scale_factor = qmmm_env%fist_scale_charge_link(ilink) + END IF + charge(i) = charge(i) * scale_factor + END IF + END IF + + CALL set_potential(potential=elp_potential,& + qeff=charge(i)) + CALL set_atomic_kind(atomic_kind=atomic_kind,& + elp_potential=elp_potential) + END DO + DEALLOCATE (charge,STAT=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CALL timestop(handle2) + !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- + ! 13. Set spline_bond_env, use for ONFO interaction !----------------------------------------------------------------------------- CALL timeset(routineN//'_onfo2',handle2) CALL pair_potential_pp_create ( potparm_bond, & @@ -2868,8 +2958,7 @@ CONTAINS atomic_kind_set , ewald_env, ff_type%ei_scale14,ff_type%vdw_scale14,& ff_type%eps_spline,max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb,& emax_spline=ff_type%emax_spline, iw=iw2, iw2=iw3, iw3=iw4,& - do_14=.TRUE., error=error) - ! point potparm ( ikind, jkind ) % spline_data => spline_data_p ( n ) + do_14=.TRUE., scale_cutoff=ff_type%scale_cutoff, error=error) DO ikind = 1, SIZE ( potparm_bond%pot, 1 ) DO jkind = ikind, SIZE ( potparm_bond%pot, 2 ) n = spline_bond_env % spltab ( ikind, jkind ) @@ -2883,7 +2972,7 @@ CONTAINS CALL timestop(handle2) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 13. ONFO + ! 14. ONFO !----------------------------------------------------------------------------- CALL timeset(routineN//'_onfo3',handle2) DO i=1,SIZE(molecule_kind_set) @@ -2923,95 +3012,6 @@ CONTAINS CALL timestop(handle2) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- - ! 14. Set up atomic_kind_set()%elp_potentail%qeff (PART 2) - !----------------------------------------------------------------------------- - CALL timeset(routineN//'_qeff',handle2) - ALLOCATE(charge(SIZE(atomic_kind_set)),STAT=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - DO i=1,SIZE(atomic_kind_set) - atomic_kind => atomic_kind_set(i) - CALL get_atomic_kind(atomic_kind=atomic_kind,& - elp_potential=elp_potential,& - atom_list=my_atom_list,& - name=atmname) - CALL get_potential(potential=elp_potential,& - qeff=charge(i)) - - found = .FALSE. - only_qm = qmmm_ff_precond_only_qm(id1=atmname,is_link=is_link_atom) - CALL uppercase(atmname) - IF(charge(i)/=-HUGE(0.0_dp)) found = .TRUE. - ! loop over params from charmm - IF(ASSOCIATED(ff_type%chm_info%charge_atm)) THEN - DO j=1,SIZE(ff_type%chm_info%charge_atm) - IF(TRIM(ff_type%chm_info%charge_atm(j))==atmname) THEN - charge(i) = ff_type%chm_info%charge(j) - IF (found) CALL Pwarning("Multiple charge declarations: "//TRIM(atmname)//& - " overwriting!",globenv,error) - found = .TRUE. - END IF - END DO - END IF - ! always have the input param last to overwrite all the other ones - IF(ASSOCIATED(ff_type%inp_info%charge_atm)) THEN - DO j=1,SIZE(ff_type%inp_info%charge_atm) - IF (iw>0) WRITE(iw,*)"Charge Checking ::",TRIM(ff_type%inp_info%charge_atm(j)),atmname - IF(TRIM(ff_type%inp_info%charge_atm(j))==atmname) THEN - charge(i) = ff_type%inp_info%charge(j) - IF (found) CALL Pwarning("Multiple charge declarations: "//TRIM(atmname)//& - " overwriting!",globenv,error) - found = .TRUE. - END IF - END DO - END IF - IF (ASSOCIATED(ff_type%inp_info%nonbonded)) THEN - IF (ASSOCIATED(ff_type%inp_info%nonbonded%pot)) THEN - ! Let's find the nonbonded potential where this atom is involved - DO j = 1, SIZE(ff_type%inp_info%nonbonded%pot) - IF (atmname==ff_type%inp_info%nonbonded%pot(j)%pot %at1.OR.& - atmname==ff_type%inp_info%nonbonded%pot(j)%pot %at2) THEN - SELECT CASE(ff_type%inp_info%nonbonded%pot(j)%pot%type) - CASE (ea_type) - ! Nothing to complain.. we don't need charges - CASE default - IF(.NOT.found) THEN - CALL stop_program("FF_coordinate_pack","missing qeff parm "//atmname) - END IF - END SELECT - END IF - END DO - END IF - END IF - ! - ! QM/MM modifications - ! - IF (only_qm.AND.my_qmmm) THEN - IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN - scale_factor = 0.0_dp - IF (is_link_atom) THEN - ! - ! Find the scaling factor... - ! - DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) - IF (ANY(my_atom_list == qmmm_env%mm_link_atoms(ilink))) EXIT - END DO - CPPostcondition(ilink <= SIZE(qmmm_env%mm_link_atoms),cp_failure_level,routineP,error,failure) - scale_factor = qmmm_env%fist_scale_charge_link(ilink) - END IF - charge(i) = charge(i) * scale_factor - END IF - END IF - - CALL set_potential(potential=elp_potential,& - qeff=charge(i)) - CALL set_atomic_kind(atomic_kind=atomic_kind,& - elp_potential=elp_potential) - END DO - DEALLOCATE (charge,STAT=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - CALL timestop(handle2) - !----------------------------------------------------------------------------- - !----------------------------------------------------------------------------- ! 15. Set spline_nonbond_env !----------------------------------------------------------------------------- CALL timeset(routineN//'_nonbond',handle2) @@ -3139,7 +3139,7 @@ CONTAINS atomic_kind_set,ewald_env,eps_spline=ff_type%eps_spline,& max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb,& emax_spline=ff_type%emax_spline, iw=iw2, iw2=iw3,iw3=iw4,& - error=error) + scale_cutoff=ff_type%scale_cutoff, error=error) !--------------------------------------- ! create the manybody spline environment !--------------------------------------- diff --git a/src/input_cp2k_mm.F b/src/input_cp2k_mm.F index 7f0a84b746..41ddc17f82 100644 --- a/src/input_cp2k_mm.F +++ b/src/input_cp2k_mm.F @@ -289,6 +289,13 @@ CONTAINS usage="EI_SCALE14 1.0", default_r_val=0.0_dp,error=error) CALL section_add_keyword(section,keyword,error=error) CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="SCALE_CUTOFF",& + description="Enables the scaling of the energy at the cutoff radius."//& + " In this way the energy at the cutoff radius is zero.",& + usage="SCALE_CUTOFF ", default_l_val=.TRUE.,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) ! ! subsections ! diff --git a/src/pair_potential.F b/src/pair_potential.F index 1b594c84ab..134d2f1b7f 100644 --- a/src/pair_potential.F +++ b/src/pair_potential.F @@ -57,6 +57,7 @@ MODULE pair_potential IMPLICIT NONE PRIVATE + REAL(KIND=dp), PARAMETER :: ifourpi = 1.0_dp / ( 4.0_dp * pi ) PUBLIC :: potential_f, potential_s PUBLIC :: spline_nonbond_control, get_nonbond_storage CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pair_potential' @@ -161,9 +162,6 @@ CONTAINS REAL(KIND=dp) :: ener_coul LOGICAL, INTENT(IN) :: my_do_14 - REAL(KIND=dp), PARAMETER :: eps = 1.0E-6_dp, & - ifourpi = 1.0_dp / ( 4.0_dp * pi ) - INTEGER :: ewald_type REAL(KIND=dp) :: alpha, eps0 @@ -179,7 +177,10 @@ CONTAINS ener_coul = ei_scale14 * ifourpi / eps0 * qi * qj / r - pot % e_cutoff_coul END IF ELSE - ener_coul = ei_scale14 * ifourpi / eps0 * qi * qj / r - pot % e_cutoff_coul + ! to avoid a large numberof points in the spline generation + ! we evaluate the 1-4 interactions analytically.. + ener_coul = 0.0_dp + !ei_scale14 * ifourpi / eps0 * qi * qj / r - pot % e_cutoff_coul ENDIF END FUNCTION ener_coul @@ -204,7 +205,7 @@ CONTAINS !****************************************************************************** SUBROUTINE spline_nonbond_control ( spline_env, potparm, atomic_kind_set, & ewald_env, ei_scale14, vdw_scale14, eps_spline, max_energy, rlow_nb,& - emax_spline, iw, iw2, iw3, do_14, error) + emax_spline, iw, iw2, iw3, do_14, scale_cutoff, error) TYPE(spline_environment_type), POINTER :: spline_env TYPE(pair_potential_pp_type), POINTER :: potparm @@ -217,6 +218,7 @@ CONTAINS emax_spline INTEGER, INTENT(IN) :: iw, iw2, iw3 LOGICAL, INTENT(IN),OPTIONAL :: do_14 + LOGICAL, INTENT(IN) :: scale_cutoff TYPE(cp_error_type), INTENT(inout), & OPTIONAL :: error @@ -226,7 +228,8 @@ CONTAINS INTEGER :: nsize, nkx, kx, ncount REAL(KIND=dp) :: dx2, e, & hicut, locut, qi, qj, x, x2,& - xdum,diffmax, xdum1, xsav, ediff, my_ei_scale14, my_vdw_scale14 + xdum,diffmax, xdum1, xsav, ediff + REAL(KIND=dp):: my_ei_scale14, my_vdw_scale14, eps0, myfac1, myfac2 TYPE(atomic_kind_type), POINTER :: atomic_kind LOGICAL :: failure, check_bounds, my_do_14 @@ -257,9 +260,11 @@ CONTAINS IF (ABS(hicut) <= 1.0E-15_dp) hicut = 1.0E3_dp potparm%pot(i,j)%pot %energy_cutoff = 0.0_dp potparm%pot(i,j)%pot %e_cutoff_coul = 0.0_dp - potparm%pot(i,j)%pot %energy_cutoff = ener_pot(potparm%pot(i,j)%pot,qi,qj,hicut,my_vdw_scale14) - potparm%pot(i,j)%pot %e_cutoff_coul & - = ener_coul ( potparm%pot( i, j )%pot, qi, qj, hicut, ewald_env, my_ei_scale14, my_do_14 ) + IF (scale_cutoff) THEN + potparm%pot(i,j)%pot %energy_cutoff = ener_pot(potparm%pot(i,j)%pot,qi,qj,hicut,my_vdw_scale14) + potparm%pot(i,j)%pot %e_cutoff_coul & + = ener_coul ( potparm%pot( i, j )%pot, qi, qj, hicut, ewald_env, my_ei_scale14, my_do_14 ) + END IF ! Find the real locut according emax_spline dx2 = (hicut-locut)/10000.0_dp x = hicut @@ -286,7 +291,18 @@ CONTAINS spline_env % spline_data_p ( n ) % spline_data % y(jx) = e x2 = x2 + dx2 END DO ! jx: loop over distance - CALL init_spline(spline_env%spline_data_p(n)%spline_data,dx=dx2) + IF (my_do_14) THEN + CALL ewald_env_get ( ewald_env, eps0 = eps0 ) + myfac1 = my_ei_scale14 * ifourpi / eps0 * qi * qj + myfac2 = 0.0_dp + IF (scale_cutoff) THEN + myfac2 = my_ei_scale14 * ifourpi / eps0 * qi * qj / hicut + END IF + CALL init_spline(spline_env%spline_data_p(n)%spline_data,dx=dx2,& + ei_scale14=myfac1, ei_cutoff14=myfac2) + ELSE + CALL init_spline(spline_env%spline_data_p(n)%spline_data,dx=dx2) + END IF ! This is the check for required accuracy on spline setup dx2 = (hicut - locut)/REAL(5*npoints+1,KIND=dp) x2 = locut + dx2 diff --git a/src/splines.F b/src/splines.F index 790d9a0e43..3ad7f3dd5f 100644 --- a/src/splines.F +++ b/src/splines.F @@ -77,6 +77,8 @@ MODULE splines ! 1/6 otherwise REAL (KIND=dp) :: x1 ! starting x value if uniform incr. REAL (KIND=dp) :: xn ! end x value if uniform incr. + ! used only for 1-4 bond_splines + REAL (KIND=dp) :: ei_scale14, ei_cutoff14 END TYPE spline_data_type TYPE spline_data_p_type @@ -378,8 +380,9 @@ CONTAINS spline_data_dest%invh = spline_data_source%invh spline_data_dest%h26 = spline_data_source%h26 spline_data_dest%x1 = spline_data_source%x1 - spline_data_dest%xn = spline_data_source%xn - + spline_data_dest%xn = spline_data_source%xn + spline_data_dest%ei_scale14 = spline_data_source%ei_scale14 + spline_data_dest%ei_cutoff14 = spline_data_source%ei_cutoff14 IF (ASSOCIATED(spline_data_source%x)) THEN ALLOCATE(spline_data_dest%x(SIZE(spline_data_source%x)), stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) @@ -535,11 +538,11 @@ CONTAINS !! SOURCE !****************************************************************************** - SUBROUTINE init_spline ( spl, dx, y1a, y1b ) + SUBROUTINE init_spline ( spl, dx, y1a, y1b, ei_scale14, ei_cutoff14) TYPE(spline_data_type), POINTER :: spl REAL(KIND=dp), INTENT(IN) :: dx - REAL(KIND=dp), INTENT(IN), OPTIONAL :: y1a, y1b + REAL(KIND=dp), INTENT(IN), OPTIONAL :: y1a, y1b, ei_scale14, ei_cutoff14 INTEGER :: err, i, n REAL(KIND=dp) :: p, s @@ -559,6 +562,10 @@ CONTAINS spl % invh = 1.0_dp/dx spl % h26 = dx**2 / 6.0_dp spl % y2 => spl % x + spl % ei_scale14 = HUGE(0.0_dp) + spl % ei_cutoff14 = HUGE(0.0_dp) + IF (PRESENT(ei_scale14 )) spl % ei_scale14 = ei_scale14 + IF (PRESENT(ei_cutoff14)) spl % ei_cutoff14 = ei_cutoff14 NULLIFY (spl%x) ALLOCATE (ww(1:n),STAT=err) IF ( err /= 0 ) CALL stop_memory ( 'init_spline', 'ww', spl % n ) diff --git a/src/topology_util.F b/src/topology_util.F index 43d1388d46..7c052805cd 100644 --- a/src/topology_util.F +++ b/src/topology_util.F @@ -2281,18 +2281,8 @@ CONTAINS atm_b = conn_info%c_bond_b(j) IF(first<=atm_a .AND. atm_a <=last) THEN DO k=1,nonfo - IF( ( (atm_a==conn_info%phi_a(k)) .AND.& - (atm_b==conn_info%phi_b(k)) ) .OR.& - ( (atm_a==conn_info%phi_b(k)) .AND.& - (atm_b==conn_info%phi_c(k)) ) .OR.& - ( (atm_a==conn_info%phi_c(k)) .AND.& - (atm_b==conn_info%phi_d(k)) ) .OR.& - ( (atm_a==conn_info%phi_d(k)) .AND.& - (atm_b==conn_info%phi_c(k)) ) .OR.& - ( (atm_a==conn_info%phi_c(k)) .AND.& - (atm_b==conn_info%phi_b(k)) ) .OR.& - ( (atm_a==conn_info%phi_b(k)) .AND.& - (atm_b==conn_info%phi_a(k)) ) ) THEN + atm_b = conn_info%onfo_a(k) + IF (atm_b >= first .AND. atm_b <= last ) THEN checkme(k)=1 END IF END DO @@ -2407,7 +2397,6 @@ CONTAINS TYPE(connectivity_info_type), POINTER :: conn_info TYPE(atom_bond_list_type), DIMENSION(:), POINTER :: ex_bond_list,& ex_bend_list,& - ex_torsion_list,& ex_onfo_list TYPE(cp_logger_type), POINTER :: logger TYPE(elp_potential_type), POINTER :: elp_potential @@ -2652,16 +2641,6 @@ CONTAINS N = 0 IF(ASSOCIATED(conn_info%theta_a)) N = SIZE(conn_info%theta_a) CALL reorder_structure(ex_bend_list, conn_info%theta_a, conn_info%theta_c, N, error) - ! Reorder torsions - ALLOCATE(ex_torsion_list(natom),stat=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - DO I=1,natom - ALLOCATE(ex_torsion_list(I)%bonds(0),stat=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - ENDDO - N = 0 - IF(ASSOCIATED(conn_info%phi_a)) N = SIZE(conn_info%phi_a) - CALL reorder_structure(ex_torsion_list, conn_info%phi_a, conn_info%phi_d, N, error) ! Reorder onfo ALLOCATE(ex_onfo_list(natom),stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) @@ -2685,7 +2664,7 @@ CONTAINS IF (topology%exclude_vdw==do_skip_13.OR. & topology%exclude_vdw==do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%bonds) dim2 = dim1 + dim2 - dim3 = SIZE(ex_torsion_list(iatom)%bonds) ! Always exclude 1-4 + dim3 = SIZE(ex_onfo_list(iatom)%bonds) ! Always exclude 1-4 dim3 = dim2 + dim3 ! They're treated in a different way.. IF (dim3 /= 0) THEN NULLIFY(list, wlist) @@ -2694,7 +2673,7 @@ CONTAINS wlist( dim0:dim0) = iatom IF (dim1>dim0) wlist(dim0+1:dim1) = ex_bond_list(iatom)%bonds IF (dim2>dim1) wlist(dim1+1:dim2) = ex_bend_list(iatom)%bonds - IF (dim3>dim2) wlist(dim2+1:dim3) = ex_torsion_list(iatom)%bonds + IF (dim3>dim2) wlist(dim2+1:dim3) = ex_onfo_list(iatom)%bonds ! Get a unique list DO i = 1, SIZE(wlist)-1 IF (wlist(i)==0)CYCLE @@ -2729,7 +2708,7 @@ CONTAINS IF (topology%exclude_ei==do_skip_13.OR. & topology%exclude_ei==do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%bonds) dim2 = dim1 + dim2 - dim3 = SIZE(ex_torsion_list(iatom)%bonds) ! Always exclude 1-4 + dim3 = SIZE(ex_onfo_list(iatom)%bonds) ! Always exclude 1-4 dim3 = dim2 + dim3 ! They're treated in a different way.. IF (dim3 /= 0) THEN ALLOCATE(wlist(dim3),stat=stat) @@ -2737,7 +2716,7 @@ CONTAINS wlist( 1: 1) = iatom wlist( 2:dim1) = ex_bond_list(iatom)%bonds wlist(dim1+1:dim2) = ex_bend_list(iatom)%bonds - wlist(dim2+1:dim3) = ex_torsion_list(iatom)%bonds + wlist(dim2+1:dim3) = ex_onfo_list(iatom)%bonds ! Get a unique list DO i = 1, SIZE(wlist)-1 IF (wlist(i)==0)CYCLE @@ -2787,13 +2766,6 @@ CONTAINS ENDDO DEALLOCATE(ex_onfo_list,stat=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - ! deallocate torsions - DO I=1,natom - DEALLOCATE(ex_torsion_list(I)%bonds,stat=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) - ENDDO - DEALLOCATE(ex_torsion_list,stat=stat) - CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) ! deallocate bends DO I=1,natom DEALLOCATE(ex_bend_list(I)%bonds,stat=stat) diff --git a/tests/Fist/regtest/TEST_FILES b/tests/Fist/regtest/TEST_FILES index fc284d1da4..b3d32cc96f 100644 --- a/tests/Fist/regtest/TEST_FILES +++ b/tests/Fist/regtest/TEST_FILES @@ -96,3 +96,5 @@ H2O-1.inp 2 # AMBER tests comparing energies w.r.t. AMBER.. differences of the order of 10^-6-10^-7 Hartree # see the h2o2_amber.inp for comments h2o2_amber.inp 2 +gly_amber.inp 2 +mol_amber.inp 2 diff --git a/tests/Fist/regtest/TEST_FILES_RESET b/tests/Fist/regtest/TEST_FILES_RESET index 8079b801e7..7965b25900 100644 --- a/tests/Fist/regtest/TEST_FILES_RESET +++ b/tests/Fist/regtest/TEST_FILES_RESET @@ -179,3 +179,7 @@ silicon_cluster_5.inp pentadiene.inp cyhex.inp cycbut.inp +#bug fix for onfo +loop.inp +#conforming the input to the amber calculation +h2o2_amber.inp diff --git a/tests/Fist/regtest/gly.psf b/tests/Fist/regtest/gly.psf new file mode 100644 index 0000000000..937bc506b9 --- /dev/null +++ b/tests/Fist/regtest/gly.psf @@ -0,0 +1,83 @@ +PSF + + 1 !NTITLE + Conversion from AMBER PARMTOP ::ACE + + 19 !NATOM + 1 MOL01 1 ACE HH31 HC 0.11230 1.00800 0 + 2 MOL01 1 ACE CH3 CT1 -0.36620 12.01000 0 + 3 MOL01 1 ACE HH32 HC 0.11230 1.00800 0 + 4 MOL01 1 ACE HH33 HC 0.11230 1.00800 0 + 5 MOL01 1 ACE C C1 0.59720 12.01000 0 + 6 MOL01 1 ACE O O1 -0.56790 16.00000 0 + 7 MOL01 2 GLY N N1 -0.37428 14.01000 0 + 8 MOL01 2 GLY H H11 0.25398 1.00800 0 + 9 MOL01 2 GLY CA CT2 -0.12884 12.01000 0 + 10 MOL01 2 GLY HA2 H0 0.08886 1.00800 0 + 11 MOL01 2 GLY HA3 H0 0.08886 1.00800 0 + 12 MOL01 2 GLY C C2 0.58058 12.01000 0 + 13 MOL01 2 GLY O O2 -0.50916 16.00000 0 + 14 MOL01 3 NME N N2 -0.41570 14.01000 0 + 15 MOL01 3 NME H H2 0.27190 1.00800 0 + 16 MOL01 3 NME CH3 CT3 -0.14900 12.01000 0 + 17 MOL01 3 NME HH31 H12 0.09760 1.00800 0 + 18 MOL01 3 NME HH32 H12 0.09760 1.00800 0 + 19 MOL01 3 NME HH33 H12 0.09760 1.00800 0 + + 18 !NBOND + 2 3 2 4 1 2 9 10 + 9 11 7 8 16 17 16 18 + 16 19 14 15 5 6 5 7 + 2 5 12 13 12 14 9 12 + 7 9 14 16 + + 30 !NTHETA + 5 7 8 4 2 5 3 2 4 + 3 2 5 1 2 3 1 2 4 + 1 2 5 12 14 15 11 9 12 + 10 9 11 10 9 12 8 7 9 + 7 9 10 7 9 11 18 16 19 + 17 16 18 17 16 19 15 14 16 + 14 16 17 14 16 18 14 16 19 + 6 5 7 5 7 9 2 5 6 + 2 5 7 13 12 14 12 14 16 + 9 12 13 9 12 14 7 9 12 + + 51 !NPHI + 6 5 7 8 6 5 7 8 + 5 7 9 10 5 7 9 10 + 5 7 9 11 5 7 9 11 + 4 2 5 6 4 2 5 6 + 4 2 5 7 3 2 5 6 + 3 2 5 6 3 2 5 6 + 3 2 5 7 2 5 7 8 + 1 2 5 6 1 2 5 6 + 1 2 5 6 1 2 5 7 + 13 12 14 15 13 12 14 15 + 12 14 16 17 12 14 16 18 + 12 14 16 19 11 9 12 13 + 11 9 12 14 11 9 12 14 + 10 9 12 13 10 9 12 14 + 10 9 12 14 9 12 14 15 + 8 7 9 10 8 7 9 11 + 8 7 9 12 15 14 16 17 + 15 14 16 18 15 14 16 19 + 5 9 7 8 12 16 14 15 + 6 5 7 9 5 7 9 12 + 5 7 9 12 5 7 9 12 + 2 5 7 9 13 12 14 16 + 9 12 14 16 7 9 12 13 + 7 9 12 14 7 9 12 14 + 7 9 12 14 2 7 5 6 + 9 14 12 13 + + 0 !NIMPHI + + 0 !NDON + + 0 !NACC + + 0 !NNB + + 0 !NGRP + diff --git a/tests/Fist/regtest/gly_amber.inp b/tests/Fist/regtest/gly_amber.inp new file mode 100644 index 0000000000..e86757932f --- /dev/null +++ b/tests/Fist/regtest/gly_amber.inp @@ -0,0 +1,61 @@ + +&FORCE_EVAL + &MM + &FORCEFIELD + ! PARM91 + ! VDW_SCALE14 0.5 + ! EI_SCALE14 0.5 + ! PARM94 and Following + VDW_SCALE14 0.5 + EI_SCALE14 0.8333333 + parmfile ../sample_pot/gly.pot + parmtype CHM + scale_cutoff .false. + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE ewald + ALPHA .36 + GMAX 51 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 30.0 30.0 30.0 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + CONN_FILE gly.psf + CONNECTIVITY PSF + COORD_FILE_NAME ../sample_pdb/gly.pdb + COORDINATE PDB + &END TOPOLOGY + &PRINT + &TOPOLOGY_INFO + UTIL_INFO + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PRINT_LEVEL DEBUG + PROGRAM FIST + PROJECT mol + RUN_TYPE md +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 1 + TIMESTEP 1.0 + TEMPERATURE 0.0 + &END MD +&END MOTION +# Amber results +# BOND = 0.0203 ANGLE = 0.0722 DIHED = 10.8805 +# VDWAALS = -0.1629 EEL = -82.1936 HBOND = 0.0000 +# 1-4 VDW = 3.5047 1-4 EEL = 48.3910 RESTRAINT = 0.0000 + + diff --git a/tests/Fist/regtest/h2o2_amber.inp b/tests/Fist/regtest/h2o2_amber.inp index e0c9ef3466..0c2b0710f7 100644 --- a/tests/Fist/regtest/h2o2_amber.inp +++ b/tests/Fist/regtest/h2o2_amber.inp @@ -10,6 +10,7 @@ EI_SCALE14 0.8333333 parmfile ../sample_pot/h2o2_amber.pot parmtype CHM + scale_cutoff .false. &END FORCEFIELD &POISSON &EWALD @@ -56,8 +57,4 @@ # BOND = 2.2271 ANGLE = 0.3320 DIHED = 0.6088 # VDWAALS = -0.0014 EEL = -0.0159 HBOND = 0.0000 # 1-4 VDW = 0.0000 1-4 EEL = 20.9514 RESTRAINT = 0.0000 -# The value of 1-4 EEL is 5 Kcal smaller because we subtract the e_cutoff_cul -# that for pure electrostatic is quite big.. contrary to AMBER -# it looks to me more nice to have the potential equals 0 at the cutoff -# AMBER does have a step at the cutoff distance diff --git a/tests/Fist/regtest/mol.psf b/tests/Fist/regtest/mol.psf new file mode 100644 index 0000000000..d4f0ec76d4 --- /dev/null +++ b/tests/Fist/regtest/mol.psf @@ -0,0 +1,36 @@ +PSF + + 1 !NTITLE + Conversion from AMBER PARMTOP ::MOL + + 6 !NATOM + 1 MOL01 1 MOL C1 C 1.00000 12.01000 0 + 2 MOL01 1 MOL O2 O -1.00000 16.00000 0 + 3 MOL01 1 MOL H3 CT 1.00000 12.01000 0 + 4 MOL01 1 MOL N4 N -3.00000 14.01000 0 + 5 MOL01 1 MOL H5 H 1.00000 1.00800 0 + 6 MOL01 1 MOL H6 H 1.00000 1.00800 0 + + 5 !NBOND + 1 3 4 5 4 6 1 2 + 1 4 + + 6 !NTHETA + 1 4 5 1 4 6 2 1 3 + 3 1 4 5 4 6 2 1 4 + + 6 !NPHI + 2 1 4 5 2 1 4 6 + 3 1 4 5 3 1 4 6 + 2 1 4 3 1 5 4 6 + + 0 !NIMPHI + + 0 !NDON + + 0 !NACC + + 0 !NNB + + 0 !NGRP + diff --git a/tests/Fist/regtest/mol_amber.inp b/tests/Fist/regtest/mol_amber.inp new file mode 100644 index 0000000000..9d02a21dd5 --- /dev/null +++ b/tests/Fist/regtest/mol_amber.inp @@ -0,0 +1,80 @@ + +&FORCE_EVAL + &MM + &FORCEFIELD + ! PARM91 + ! VDW_SCALE14 0.5 + ! EI_SCALE14 0.5 + ! PARM94 and Following + VDW_SCALE14 0.5 + EI_SCALE14 0.83333333333333333333 + parmfile ../sample_pot/mol.pot + parmtype CHM + scale_cutoff .false. +#&charge +#atom C +#charge 0.0 +#&end +#&charge +#atom O +#charge 0.0 +#&end +#&charge +#atom N +#charge 0.0 +#&end +#&charge +#atom H1 +#charge 0.0 +#&end +#&charge +#atom CT +#charge 0.0 +#&end + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE ewald + ALPHA .36 + GMAX 51 + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + ABC 30.0 30.0 30.0 + UNIT ANGSTROM + &END CELL + &TOPOLOGY + CONN_FILE mol.psf + CONNECTIVITY PSF + COORD_FILE_NAME ../sample_pdb/mol.pdb + COORDINATE PDB + &END TOPOLOGY + &PRINT + &TOPOLOGY_INFO + UTIL_INFO + &END + &END + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + FFTLIB FFTSG + PRINT_LEVEL DEBUG + PROGRAM FIST + PROJECT mol + RUN_TYPE md +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 1 + TIMESTEP 1.0 + TEMPERATURE 0.0 + &END MD +&END MOTION +# BOND = 0.0001 ANGLE = 0.0001 DIHED = 4.0000 +# VDWAALS = -0.0016 EEL = -0.1492 HBOND = 0.0000 +# 1-4 VDW = -0.0565 1-4 EEL = -8.1447 RESTRAINT = 0.0000 + + diff --git a/tests/Fist/sample_pdb/gly.pdb b/tests/Fist/sample_pdb/gly.pdb new file mode 100644 index 0000000000..4816872767 --- /dev/null +++ b/tests/Fist/sample_pdb/gly.pdb @@ -0,0 +1,22 @@ +REMARK ACE +ATOM 1 1HH3 ACE 1 2.000 1.000 0.000 +ATOM 2 CH3 ACE 1 2.000 2.090 0.000 +ATOM 3 2HH3 ACE 1 1.486 2.454 0.890 +ATOM 4 3HH3 ACE 1 1.486 2.454 -0.890 +ATOM 5 C ACE 1 3.427 2.641 0.000 +ATOM 6 O ACE 1 4.391 1.877 0.000 +ATOM 7 N GLY 2 3.555 3.970 0.000 +ATOM 8 H GLY 2 2.733 4.556 0.000 +ATOM 9 CA GLY 2 4.853 4.614 0.000 +ATOM 10 2HA GLY 2 5.408 4.316 0.890 +ATOM 11 3HA GLY 2 5.408 4.316 -0.890 +ATOM 12 C GLY 2 4.694 6.128 0.000 +ATOM 13 O GLY 2 3.576 6.637 0.000 +ATOM 14 N NME 3 5.819 6.847 0.000 +ATOM 15 H NME 3 6.716 6.382 0.000 +ATOM 16 CH3 NME 3 5.801 8.296 0.000 +ATOM 17 1HH3 NME 3 4.769 8.647 0.000 +ATOM 18 2HH3 NME 3 6.310 8.666 0.890 +ATOM 19 3HH3 NME 3 6.310 8.666 -0.890 +TER +END diff --git a/tests/Fist/sample_pdb/mol.pdb b/tests/Fist/sample_pdb/mol.pdb new file mode 100644 index 0000000000..1ae83367be --- /dev/null +++ b/tests/Fist/sample_pdb/mol.pdb @@ -0,0 +1,9 @@ +REMARK MOL +ATOM 1 C1 MOL 1 2.155 -0.393 0.000 +ATOM 2 O2 MOL 1 2.881 -1.384 0.000 +ATOM 3 H3 MOL 1 0.640 -0.546 0.000 +ATOM 4 N4 MOL 1 2.630 0.855 0.000 +ATOM 5 H5 MOL 1 1.993 1.638 0.000 +ATOM 6 H6 MOL 1 3.628 1.015 0.000 +TER +END diff --git a/tests/Fist/sample_pot/gly.pot b/tests/Fist/sample_pot/gly.pot new file mode 100644 index 0000000000..2358927b69 --- /dev/null +++ b/tests/Fist/sample_pot/gly.pot @@ -0,0 +1,165 @@ +*>>>>>>> AMBER FF Converted into CHARMM FF style <<<<<<< +*>>>>>>> Generated on :: 20060602 210843.445 +0200 by :: teo <<<<<<< +*>>>>>>> pcihpc07 <<<<<<< +*>>>>>>> Leap Title :: ACE <<<<<<< +*>>>>>>> Send all comments related to the FFs conversion to <<<<<<< +*>>>>>>> teodoro.laino@gmail.com <<<<<<< + +BONDS +! +!V(bond) = Kb(b - b0)**2 +! +!Kb: kcal/mole/A**2 +!b0: A +! +!atom type Kb b0 +! +CT1 HC 340.000000000 1.090000000 +CT2 H0 340.000000000 1.090000000 +N1 H11 434.000000000 1.010000000 +CT3 H12 340.000000000 1.090000000 +N2 H2 434.000000000 1.010000000 +C1 O1 570.000000000 1.229000000 +C1 N1 490.000000000 1.335000000 +CT1 C1 317.000000000 1.522000000 +C2 O2 570.000000000 1.229000000 +C2 N2 490.000000000 1.335000000 +CT2 C2 317.000000000 1.522000000 +N1 CT2 337.000000000 1.449000000 +N2 CT3 337.000000000 1.449000000 + +ANGLES +! +!V(angle) = Ktheta(Theta - Theta0)**2 +! +!V(Urey-Bradley) = Kub(S - S0)**2 +! +!Ktheta: kcal/mole/rad**2 +!Theta0: degrees +!Kub: kcal/mole/A**2 (Urey-Bradley) +!S0: A +! +!atom types Ktheta Theta0 Kub S0 +! +C1 N1 H11 50.000000000 120.000051429 +HC CT1 C1 50.000000000 109.500046929 +HC CT1 HC 35.000000000 109.500046929 +C2 N2 H2 50.000000000 120.000051429 +H0 CT2 C2 50.000000000 109.500046929 +H0 CT2 H0 35.000000000 109.500046929 +H11 N1 CT2 50.000000000 118.040050474 +N1 CT2 H0 50.000000000 109.500046929 +H12 CT3 H12 35.000000000 109.500046929 +H2 N2 CT3 50.000000000 118.040050474 +N2 CT3 H12 50.000000000 109.500046929 +O1 C1 N1 80.000000000 122.900052672 +C1 N1 CT2 50.000000000 121.900052243 +CT1 C1 O1 80.000000000 120.400051601 +CT1 C1 N1 70.000000000 116.600049972 +O2 C2 N2 80.000000000 122.900052672 +C2 N2 CT3 50.000000000 121.900052243 +CT2 C2 O2 80.000000000 120.400051601 +CT2 C2 N2 70.000000000 116.600049972 +N1 CT2 C2 63.000000000 110.100047186 + +DIHEDRALS +! +!V(dihedral) = Kchi(1 + cos(n(chi) - delta)) +! +!Kchi: kcal/mole +!n: multiplicity +!delta: degrees +! +!atom types Kchi n delta +! +O1 C1 N1 H11 2.000000000 1 0.000000000 +O1 C1 N1 H11 2.500000000 2 180.000077144 +C1 N1 CT2 H0 0.457500000 1 0.000000000 +C1 N1 CT2 H0 1.255800000 2 180.000077144 +HC CT1 C1 O1 0.800000000 1 0.000000000 +HC CT1 C1 O1 0.080000000 3 180.000077144 +HC CT1 C1 N1 0.000000000 2 0.000000000 +HC CT1 C1 O1 0.000000000 2 0.000000000 +CT1 C1 N1 H11 2.500000000 2 180.000077144 +O2 C2 N2 H2 2.000000000 1 0.000000000 +O2 C2 N2 H2 2.500000000 2 180.000077144 +C2 N2 CT3 H12 0.000000000 2 0.000000000 +H0 CT2 C2 O2 0.000000000 2 0.000000000 +H0 CT2 C2 N2 1.060700000 1 180.000077144 +H0 CT2 C2 N2 0.011000000 2 0.000000000 +CT2 C2 N2 H2 2.500000000 2 180.000077144 +H11 N1 CT2 H0 0.000000000 2 0.000000000 +H11 N1 CT2 C2 0.000000000 2 0.000000000 +H2 N2 CT3 H12 0.000000000 2 0.000000000 +O1 C1 N1 CT2 2.500000000 2 180.000077144 +C1 N1 CT2 C2 1.015900000 1 0.000000000 +C1 N1 CT2 C2 0.345100000 2 180.000077144 +C1 N1 CT2 C2 0.225900000 3 0.000000000 +CT1 C1 N1 CT2 2.500000000 2 180.000077144 +O2 C2 N2 CT3 2.500000000 2 180.000077144 +CT2 C2 N2 CT3 2.500000000 2 180.000077144 +N1 CT2 C2 O2 0.000000000 2 0.000000000 +N1 CT2 C2 N2 0.683900000 1 180.000077144 +N1 CT2 C2 N2 1.453700000 2 180.000077144 +N1 CT2 C2 N2 0.461500000 3 180.000077144 +C1 CT2 N1 H11 1.100000000 2 180.000077144 +C2 CT3 N2 H2 1.100000000 2 180.000077144 +CT1 N1 C1 O1 10.500000000 2 180.000077144 +CT2 N2 C2 O2 10.500000000 2 180.000077144 + +IMPROPER +! +!V(improper) = Kpsi(psi - psi0)**2 +! +!Kpsi: kcal/mole/rad**2 +!psi0: degrees +!note that the second column of numbers (0) is ignored +! +!atom types Kpsi psi0 +! + +NONBONDED +! +!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6] +! +!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j) +!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j +! +!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4 +! +HC 0.000000000 0.015700000 1.487000000 +CT1 0.000000000 0.109400000 1.908000000 +C1 0.000000000 0.086000000 1.908000000 +O1 0.000000000 0.210000000 1.661200000 +N1 0.000000000 0.170000000 1.824000000 +H11 0.000000000 0.015700000 0.600000000 +CT2 0.000000000 0.109400000 1.908000000 +H0 0.000000000 0.015700000 1.386999999 +C2 0.000000000 0.086000000 1.908000000 +O2 0.000000000 0.210000000 1.661200000 +N2 0.000000000 0.170000000 1.824000000 +H2 0.000000000 0.015700000 0.600000000 +CT3 0.000000000 0.109400000 1.908000000 +H12 0.000000000 0.015700000 1.386999999 + +END + +! +! This Section can be cutted & pasted into the Fist input file.. +! +CHARGES +HC 0.112300000 +CT1 -0.366200000 +C1 0.597200002 +O1 -0.567900002 +N1 -0.374282000 +H11 0.253981000 +CT2 -0.128844000 +H0 0.088859000 +C2 0.580583999 +O2 -0.509157000 +N2 -0.415700000 +H2 0.271900000 +CT3 -0.149000000 +H12 0.097600000 +END CHARGES diff --git a/tests/Fist/sample_pot/mol.pot b/tests/Fist/sample_pot/mol.pot new file mode 100644 index 0000000000..db53518573 --- /dev/null +++ b/tests/Fist/sample_pot/mol.pot @@ -0,0 +1,94 @@ +*>>>>>>> AMBER FF Converted into CHARMM FF style <<<<<<< +*>>>>>>> Generated on :: 20060602 210908.077 +0200 by :: teo <<<<<<< +*>>>>>>> pcihpc07 <<<<<<< +*>>>>>>> Leap Title :: MOL <<<<<<< +*>>>>>>> Send all comments related to the FFs conversion to <<<<<<< +*>>>>>>> teodoro.laino@gmail.com <<<<<<< + +BONDS +! +!V(bond) = Kb(b - b0)**2 +! +!Kb: kcal/mole/A**2 +!b0: A +! +!atom type Kb b0 +! +C CT 317.000000000 1.522000000 +N H 434.000000000 1.010000000 +C O 570.000000000 1.229000000 +C N 490.000000000 1.335000000 + +ANGLES +! +!V(angle) = Ktheta(Theta - Theta0)**2 +! +!V(Urey-Bradley) = Kub(S - S0)**2 +! +!Ktheta: kcal/mole/rad**2 +!Theta0: degrees +!Kub: kcal/mole/A**2 (Urey-Bradley) +!S0: A +! +!atom types Ktheta Theta0 Kub S0 +! +C N H 50.000000000 120.000051429 +O C CT 80.000000000 120.400051601 +CT C N 70.000000000 116.600049972 +H N H 35.000000000 120.000051429 +O C N 80.000000000 122.900052672 + +DIHEDRALS +! +!V(dihedral) = Kchi(1 + cos(n(chi) - delta)) +! +!Kchi: kcal/mole +!n: multiplicity +!delta: degrees +! +!atom types Kchi n delta +! +O C N H 2.000000000 1 0.000000000 +O C N H 2.500000000 2 180.000077144 +CT C N H 2.500000000 2 180.000077144 +O C N CT 10.500000000 2 180.000077144 +C H N H 1.000000000 2 180.000077144 + +IMPROPER +! +!V(improper) = Kpsi(psi - psi0)**2 +! +!Kpsi: kcal/mole/rad**2 +!psi0: degrees +!note that the second column of numbers (0) is ignored +! +!atom types Kpsi psi0 +! + +NONBONDED +! +!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6] +! +!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j) +!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j +! +!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4 +! +C 0.000000000 0.086000000 1.908000000 +O 0.000000000 0.210000000 1.661200000 +CT 0.000000000 0.109400000 1.908000000 +N 0.000000000 0.170000000 1.824000000 +H 0.000000000 0.015700000 0.600000000 + +END + +! +! This Section can be cutted & pasted into the Fist input file.. +! +CHARGES +C 1.000000000 +O -1.000000000 +CT 1.000000000 +N -3.000000000 +H 1.000000000 +END CHARGES diff --git a/tools/Fist/topology_converter/leap2fist.F b/tools/Fist/topology_converter/leap2fist.F index dd4fbb468a..dfcf515413 100644 --- a/tools/Fist/topology_converter/leap2fist.F +++ b/tools/Fist/topology_converter/leap2fist.F @@ -821,7 +821,7 @@ CONTAINS IMPLICIT NONE INTEGER, INTENT(IN) :: unit INTEGER :: my_sec, my_res, my_bonds, local, my_phi, my_theta, my_atom - INTEGER :: i, ind, ifirst, ilast + INTEGER :: i, j, ind, ifirst, ilast INTEGER, DIMENSION(:), POINTER :: my_i, my_j, my_k, my_l CHARACTER (LEN=2) :: int_format @@ -894,6 +894,7 @@ CONTAINS my_phi = COUNT(MASK=LPH>=0) + COUNT(MASK=LP>=0) ALLOCATE(my_i(my_phi), my_j(my_phi), my_k(my_phi), my_l(my_phi)) ind = 0 + ! Get dihedrals DO I = 1, NPHIH IF (LPH(i) >= 0) THEN ind = ind + 1 @@ -913,8 +914,41 @@ CONTAINS END IF END DO IF (ind .NE. my_phi) CALL stop_converter("Error in write_psf :: evaluation of proper torsion.") - WRITE(unit,*) " ",my_phi," !NPHI" - WRITE(unit,'(8'//int_format//')') (my_i(local),my_j(local),my_k(local),my_l(local), local=1,my_phi) + ! Now get unique dihedrals + DO i = 1, my_phi + IF ( my_i(i)==-999999 .AND.& + my_j(i)==-999999 .AND.& + my_k(i)==-999999 .AND.& + my_l(i)==-999999) CYCLE + DO j = i+1, my_phi + IF ( my_i(i)==my_i(j) .AND.& + my_j(i)==my_j(j) .AND.& + my_k(i)==my_k(j) .AND.& + my_l(i)==my_l(j)) THEN + my_i(j)=-999999 + my_j(j)=-999999 + my_k(j)=-999999 + my_l(j)=-999999 + ind = ind - 1 + END IF + END DO + END DO + DO WHILE (ANY(my_i(1:ind)==-999999)) + DO i = 1, my_phi + IF ( my_i(i)==-999999 .AND.& + my_j(i)==-999999 .AND.& + my_k(i)==-999999 .AND.& + my_l(i)==-999999) THEN + my_i(i:i+my_phi-(i+1)) = my_i(i+1:my_phi) + my_j(i:i+my_phi-(i+1)) = my_j(i+1:my_phi) + my_k(i:i+my_phi-(i+1)) = my_k(i+1:my_phi) + my_l(i:i+my_phi-(i+1)) = my_l(i+1:my_phi) + END IF + END DO + END DO + ! + WRITE(unit,*) " ",ind," !NPHI" + WRITE(unit,'(8'//int_format//')') (my_i(local),my_j(local),my_k(local),my_l(local), local=1,ind) WRITE(unit,*) "" DEALLOCATE(my_i, my_j, my_k, my_l) ! @@ -942,8 +976,40 @@ CONTAINS END IF END DO IF (ind .NE. my_phi) CALL stop_converter("Error in write_psf :: evaluation of improper torsion.") - WRITE(unit,*) " ",my_phi," !NIMPHI" - WRITE(unit,'(8'//int_format//')') (my_i(local),my_j(local),my_k(local),my_l(local), local=1,my_phi) + ! Now get unique dihedrals + DO i = 1, my_phi + IF ( my_i(i)==-999999 .AND.& + my_j(i)==-999999 .AND.& + my_k(i)==-999999 .AND.& + my_l(i)==-999999) CYCLE + DO j = i+1, my_phi + IF ( my_i(i)==my_i(j) .AND.& + my_j(i)==my_j(j) .AND.& + my_k(i)==my_k(j) .AND.& + my_l(i)==my_l(j)) THEN + my_i(j)=-999999 + my_j(j)=-999999 + my_k(j)=-999999 + my_l(j)=-999999 + ind = ind - 1 + END IF + END DO + END DO + DO WHILE (ANY(my_i(1:ind)==-999999)) + DO i = 1, my_phi + IF ( my_i(i)==-999999 .AND.& + my_j(i)==-999999 .AND.& + my_k(i)==-999999 .AND.& + my_l(i)==-999999) THEN + my_i(i:i+my_phi-(i+1)) = my_i(i+1:my_phi) + my_j(i:i+my_phi-(i+1)) = my_j(i+1:my_phi) + my_k(i:i+my_phi-(i+1)) = my_k(i+1:my_phi) + my_l(i:i+my_phi-(i+1)) = my_l(i+1:my_phi) + END IF + END DO + END DO + WRITE(unit,*) " ",ind," !NIMPHI" + WRITE(unit,'(8'//int_format//')') (my_i(local),my_j(local),my_k(local),my_l(local), local=1,ind) WRITE(unit,*) "" DEALLOCATE(my_i, my_j, my_k, my_l) ELSE @@ -965,8 +1031,40 @@ CONTAINS my_l(ind) = ABS(LP(i))/3+1 END DO IF (ind .NE. my_phi) CALL stop_converter("Error in write_psf :: evaluation of proper torsion.") - WRITE(unit,*) " ",my_phi," !NPHI" - WRITE(unit,'(8'//int_format//')') (my_i(local),my_j(local),my_k(local),my_l(local), local=1,my_phi) + ! Now get unique dihedrals + DO i = 1, my_phi + IF ( my_i(i)==-999999 .AND.& + my_j(i)==-999999 .AND.& + my_k(i)==-999999 .AND.& + my_l(i)==-999999) CYCLE + DO j = i+1, my_phi + IF ( my_i(i)==my_i(j) .AND.& + my_j(i)==my_j(j) .AND.& + my_k(i)==my_k(j) .AND.& + my_l(i)==my_l(j)) THEN + my_i(j)=-999999 + my_j(j)=-999999 + my_k(j)=-999999 + my_l(j)=-999999 + ind = ind - 1 + END IF + END DO + END DO + DO WHILE (ANY(my_i(1:ind)==-999999)) + DO i = 1, my_phi + IF ( my_i(i)==-999999 .AND.& + my_j(i)==-999999 .AND.& + my_k(i)==-999999 .AND.& + my_l(i)==-999999) THEN + my_i(i:i+my_phi-(i+1)) = my_i(i+1:my_phi) + my_j(i:i+my_phi-(i+1)) = my_j(i+1:my_phi) + my_k(i:i+my_phi-(i+1)) = my_k(i+1:my_phi) + my_l(i:i+my_phi-(i+1)) = my_l(i+1:my_phi) + END IF + END DO + END DO + WRITE(unit,*) " ",ind," !NPHI" + WRITE(unit,'(8'//int_format//')') (my_i(local),my_j(local),my_k(local),my_l(local), local=1,ind) WRITE(unit,*) "" DEALLOCATE(my_i, my_j, my_k, my_l) WRITE(unit,*) " ",0," !NIMPHI"