comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-26 13:15:21 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Bug fix for onfo. Fully conforming amber results to FIST results. introduced few

Browse source 
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
This commit is contained in:
Teodoro Laino 2006-06-02 17:27:06 +00:00
parent f9fef3c64b
commit 71b2f738d4
19 changed files with 815 additions and 154 deletions
Download patch file Download diff file Expand all files Collapse all files

9
src/fist_intra_force.F
View file

@ -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

1
src/force_field_types.F
View file

@ -168,6 +168,7 @@ MODULE force_field_types
max_energy,&
emax_spline,&
rlow_nb
LOGICAL :: scale_cutoff
END TYPE force_field_type
! *** Public subroutines ***

188
src/force_fields.F
View file

@ -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
!---------------------------------------

7
src/input_cp2k_mm.F
View file

@ -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 <LOGICAL>", default_l_val=.TRUE.,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
!
! subsections
!

36
src/pair_potential.F
View file

@ -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

15
src/splines.F
View file

@ -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 )

40
src/topology_util.F
View file

@ -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)

2
tests/Fist/regtest/TEST_FILES
View file

@ -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

4
tests/Fist/regtest/TEST_FILES_RESET
View file

@ -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

83
tests/Fist/regtest/gly.psf Normal file
View file

@ -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

61
tests/Fist/regtest/gly_amber.inp Normal file
View file

@ -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

5
tests/Fist/regtest/h2o2_amber.inp
View file

@ -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

36
tests/Fist/regtest/mol.psf Normal file
View file

@ -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

80
tests/Fist/regtest/mol_amber.inp Normal file
View file

@ -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

22
tests/Fist/sample_pdb/gly.pdb Normal file
View file

@ -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

9
tests/Fist/sample_pdb/mol.pdb Normal file
View file

@ -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

165
tests/Fist/sample_pot/gly.pot Normal file
View file

@ -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

94
tests/Fist/sample_pot/mol.pot Normal file
View file

@ -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

112
tools/Fist/topology_converter/leap2fist.F
View file

@ -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"