comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

correct the analytical form of the dihedral potential

Browse source 
svn-origin-rev: 7903
This commit is contained in:
Teodoro Laino 2008-11-07 18:09:18 +00:00
parent 64521c5416
commit 6aac6bb711
24 changed files with 178 additions and 151 deletions
Download patch file Download diff file Expand all files Collapse all files

42
src/fist_intra_force.F
View file

@ -85,12 +85,11 @@ CONTAINS
index_a, index_b, index_c, index_d, ionef, ishell, itorsion, nbends, &
nbonds, nimptors, nkind, nmol_per_kind, nonefs, nshell, ntorsions, nub
REAL(KIND=dp) :: d12, d32, dist, dist1, dist2, energy, fscalar, id12, &
id32, is32, ism, isn, iw1, iw2, rab2, rinv, rinv2, s32, sm, sn, theta, &
w1, w2
id32, is32, ism, isn, rab2, rinv, rinv2, s32, sm, sn, theta
REAL(KIND=dp), DIMENSION(3) :: b12, b32, g1, g2, g3, gt1, &
gt2, gt3, gt4, k1, k2, k3, &
k4, rij, t12, t32, t34, t43, &
tm, tn, u1, u2
tm, tn
TYPE(bend_type), POINTER :: bend_list( : )
TYPE(bond_type), POINTER :: bond_list( : )
TYPE(cp_logger_type), POINTER :: logger
@ -310,26 +309,33 @@ CONTAINS
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
t34 = particle_set ( index_c ) % r - particle_set ( index_d ) % r
t43 = particle_set ( index_d ) % r - particle_set ( index_c ) % r
t12 = pbc(t12,cell)
t32 = pbc(t32,cell)
t34 = pbc(t34,cell)
t43 = pbc(t43,cell)
s32 = SQRT(DOT_PRODUCT(t32,t32))
! t12 x t32
tm(1) = t12(2)*t32(3)-t32(2)*t12(3)
tm(2) = -t12(1)*t32(3)+t32(1)*t12(3)
tm(3) = t12(1)*t32(2)-t32(1)*t12(2)
! t32 x t34
tn(1) = t32(2)*t34(3)-t34(2)*t32(3)
tn(2) = -t32(1)*t34(3)+t34(1)*t32(3)
tn(3) = t32(1)*t34(2)-t34(1)*t32(2)
sm = SQRT ( DOT_PRODUCT(tm,tm) )
ism = 1.0_dp/sm
sn = SQRT ( DOT_PRODUCT(tn,tn) )
isn = 1.0_dp/sn
s32 = SQRT ( DOT_PRODUCT(t32,t32) )
is32 = 1.0_dp/s32
dist1 = DOT_PRODUCT(t12,t32)
dist2 = DOT_PRODUCT(t43,t32)
u1 = t12 - dist1*t32*is32**2
u2 = t43 - dist2*t32*is32**2
w1 = SQRT(DOT_PRODUCT(u1,u1))
w2 = SQRT(DOT_PRODUCT(u2,u2))
IF ((ABS(w1)<EPSILON(0.0_dp)).OR.(ABS(w2)<EPSILON(0.0_dp))) CYCLE
iw1 = 1.0_dp/w1
iw2 = 1.0_dp/w2
dist2 = DOT_PRODUCT(t34,t32)
DO imul = 1, torsion_list(itorsion)%torsion_kind%nmul
CALL force_torsions(torsion_list(itorsion)%torsion_kind%id_type, &
is32,dist1,dist2,u1,u2,iw1,iw2,&
s32, is32,ism,isn,dist1,dist2,tm,tn, t12, t32,&
torsion_list(itorsion)%torsion_kind%k(imul), &
torsion_list(itorsion)%torsion_kind%cosphi0(imul), &
torsion_list(itorsion)%torsion_kind%phi0(imul), &
torsion_list(itorsion)%torsion_kind%m(imul), &
gt1,gt2,gt3,gt4,energy,fscalar)
pot_torsion = pot_torsion + energy
@ -376,13 +382,13 @@ CONTAINS
t34 = pbc(t34,cell)
t43 = pbc(t43,cell)
! t12 x t32
tm(1) = t12(2)*t32(3)-t32(2)*t12(3)
tm(1) = t12(2)*t32(3)-t32(2)*t12(3)
tm(2) = -t12(1)*t32(3)+t32(1)*t12(3)
tm(3) = t12(1)*t32(2)-t32(1)*t12(2)
tm(3) = t12(1)*t32(2)-t32(1)*t12(2)
! t32 x t34
tn(1) = t32(2)*t34(3)-t34(2)*t32(3)
tn(1) = t32(2)*t34(3)-t34(2)*t32(3)
tn(2) = -t32(1)*t34(3)+t34(1)*t32(3)
tn(3) = t32(1)*t34(2)-t34(1)*t32(2)
tn(3) = t32(1)*t34(2)-t34(1)*t32(2)
sm = SQRT ( DOT_PRODUCT(tm,tm) )
ism = 1.0_dp/sm
sn = SQRT ( DOT_PRODUCT(tn,tn) )

40
src/force_field_types.F
View file

@ -62,7 +62,7 @@ MODULE force_field_types
INTEGER :: id_type
INTEGER :: nmul
INTEGER, POINTER :: m(:)
REAL(KIND = dp), POINTER :: k(:),cosphi0(:)
REAL(KIND = dp), POINTER :: k(:),phi0(:)
INTEGER :: kind_number
END TYPE torsion_kind_type
@ -115,7 +115,7 @@ MODULE force_field_types
CHARACTER (LEN=default_string_length), POINTER :: torsion_d(:)
REAL(KIND=dp), POINTER :: torsion_k(:)
INTEGER, POINTER :: torsion_m(:)
REAL(KIND=dp), POINTER :: torsion_cosphi0(:)
REAL(KIND=dp), POINTER :: torsion_phi0(:)
INTEGER, POINTER :: impr_kind(:)
CHARACTER (LEN=default_string_length), POINTER :: impr_a(:)
CHARACTER (LEN=default_string_length), POINTER :: impr_b(:)
@ -150,7 +150,7 @@ MODULE force_field_types
CHARACTER (LEN=default_string_length), POINTER :: torsion_d(:)
REAL(KIND=dp), POINTER :: torsion_k(:)
INTEGER, POINTER :: torsion_m(:)
REAL(KIND=dp), POINTER :: torsion_cosphi0(:)
REAL(KIND=dp), POINTER :: torsion_phi0(:)
CHARACTER (LEN=default_string_length), POINTER :: impr_a(:)
CHARACTER (LEN=default_string_length), POINTER :: impr_b(:)
CHARACTER (LEN=default_string_length), POINTER :: impr_c(:)
@ -182,7 +182,7 @@ MODULE force_field_types
CHARACTER (LEN=default_string_length), POINTER :: torsion_d(:)
REAL(KIND=dp), POINTER :: torsion_k(:)
REAL(KIND=dp), POINTER :: torsion_m(:)
REAL(KIND=dp), POINTER :: torsion_cosphi0(:)
REAL(KIND=dp), POINTER :: torsion_phi0(:)
CHARACTER (LEN=default_string_length), POINTER :: nonbond_a(:)
REAL (KIND=dp), POINTER :: nonbond_eps(:)
REAL (KIND=dp), POINTER :: nonbond_rmin2(:)
@ -199,7 +199,7 @@ MODULE force_field_types
REAL(KIND=dp), POINTER :: bend_theta0(:)
REAL(KIND=dp), POINTER :: torsion_k(:)
INTEGER, POINTER :: torsion_m(:)
REAL(KIND=dp), POINTER :: torsion_cosphi0(:)
REAL(KIND=dp), POINTER :: torsion_phi0(:)
REAL(KIND=dp), POINTER :: impr_k(:)
REAL(KIND=dp), POINTER :: impr_phi0(:)
CHARACTER (LEN=default_string_length), POINTER :: nonbond_a(:)
@ -361,7 +361,7 @@ CONTAINS
torsion_kind_set(ikind)%nmul = 0
NULLIFY(torsion_kind_set(ikind)%k)
NULLIFY(torsion_kind_set(ikind)%m)
NULLIFY(torsion_kind_set(ikind)%cosphi0)
NULLIFY(torsion_kind_set(ikind)%phi0)
torsion_kind_set(ikind)%kind_number = ikind
END DO
END IF
@ -516,7 +516,7 @@ CONTAINS
NULLIFY(gro_info%bend_theta0)
NULLIFY(gro_info%torsion_k)
NULLIFY(gro_info%torsion_m)
NULLIFY(gro_info%torsion_cosphi0)
NULLIFY(gro_info%torsion_phi0)
NULLIFY(gro_info%impr_k)
NULLIFY(gro_info%impr_phi0)
NULLIFY(gro_info%nonbond_a)
@ -557,7 +557,7 @@ CONTAINS
NULLIFY(chm_info%torsion_d)
NULLIFY(chm_info%torsion_k)
NULLIFY(chm_info%torsion_m)
NULLIFY(chm_info%torsion_cosphi0)
NULLIFY(chm_info%torsion_phi0)
NULLIFY(chm_info%impr_a)
NULLIFY(chm_info%impr_b)
NULLIFY(chm_info%impr_c)
@ -597,7 +597,7 @@ CONTAINS
NULLIFY(amb_info%torsion_d)
NULLIFY(amb_info%torsion_k)
NULLIFY(amb_info%torsion_m)
NULLIFY(amb_info%torsion_cosphi0)
NULLIFY(amb_info%torsion_phi0)
NULLIFY(amb_info%nonbond_a)
NULLIFY(amb_info%nonbond_eps)
NULLIFY(amb_info%nonbond_rmin2)
@ -646,7 +646,7 @@ CONTAINS
NULLIFY(inp_info%torsion_d)
NULLIFY(inp_info%torsion_k)
NULLIFY(inp_info%torsion_m)
NULLIFY(inp_info%torsion_cosphi0)
NULLIFY(inp_info%torsion_phi0)
NULLIFY(inp_info%impr_kind)
NULLIFY(inp_info%impr_a)
NULLIFY(inp_info%impr_b)
@ -764,8 +764,8 @@ CONTAINS
DEALLOCATE(gro_info%torsion_m,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(gro_info%torsion_cosphi0)) THEN
DEALLOCATE(gro_info%torsion_cosphi0,STAT=stat)
IF(ASSOCIATED(gro_info%torsion_phi0)) THEN
DEALLOCATE(gro_info%torsion_phi0,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(gro_info%impr_k)) THEN
@ -900,8 +900,8 @@ CONTAINS
DEALLOCATE(chm_info%torsion_m,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(chm_info%torsion_cosphi0)) THEN
DEALLOCATE(chm_info%torsion_cosphi0,STAT=stat)
IF(ASSOCIATED(chm_info%torsion_phi0)) THEN
DEALLOCATE(chm_info%torsion_phi0,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(chm_info%impr_a)) THEN
@ -1032,8 +1032,8 @@ CONTAINS
DEALLOCATE(amb_info%torsion_m,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(amb_info%torsion_cosphi0)) THEN
DEALLOCATE(amb_info%torsion_cosphi0,STAT=stat)
IF(ASSOCIATED(amb_info%torsion_phi0)) THEN
DEALLOCATE(amb_info%torsion_phi0,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(amb_info%nonbond_a)) THEN
@ -1200,8 +1200,8 @@ CONTAINS
DEALLOCATE(inp_info%torsion_m,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(inp_info%torsion_cosphi0)) THEN
DEALLOCATE(inp_info%torsion_cosphi0,STAT=stat)
IF(ASSOCIATED(inp_info%torsion_phi0)) THEN
DEALLOCATE(inp_info%torsion_phi0,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(inp_info%impr_kind)) THEN
@ -1309,8 +1309,8 @@ CONTAINS
DEALLOCATE(torsion_kind%m,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
IF(ASSOCIATED(torsion_kind%cosphi0)) THEN
DEALLOCATE(torsion_kind%cosphi0,STAT=stat)
IF(ASSOCIATED(torsion_kind%phi0)) THEN
DEALLOCATE(torsion_kind%phi0,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END IF

2
src/force_fields.F
View file

@ -244,7 +244,7 @@ CONTAINS
DO i = 1, SIZE(ff_type%amb_info%torsion_a)
k = cp_unit_from_cp2k(ff_type%amb_info%torsion_k(i),"kcalmol",error=error)
m = ff_type%amb_info%torsion_m(i)
phi0 = cp_unit_from_cp2k(ff_type%amb_info%torsion_cosphi0(i),"deg",error=error)
phi0 = cp_unit_from_cp2k(ff_type%amb_info%torsion_phi0(i),"deg",error=error)
WRITE(iw, 2003) ff_type%amb_info%torsion_a(i),&
ff_type%amb_info%torsion_b(i),&
ff_type%amb_info%torsion_c(i),&

28
src/force_fields_all.F
View file

@ -1309,19 +1309,19 @@ CONTAINS
IF(itype>0) THEN
CALL reallocate(torsion_list(j)%torsion_kind%k,1,1)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,1)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,1)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,1)
torsion_list(j)%torsion_kind%nmul = 1
torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype)
torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype)
torsion_list(j)%torsion_kind%cosphi0(1) = gro_info%torsion_cosphi0(itype)
torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype)
ELSE
CALL reallocate(torsion_list(j)%torsion_kind%k,1,1)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,1)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,1)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,1)
torsion_list(j)%torsion_kind%nmul = 1
torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype/k)
torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype/k)
torsion_list(j)%torsion_kind%cosphi0(1) = gro_info%torsion_cosphi0(itype/k)
torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype/k)
END IF
torsion_list(j)%torsion_kind%id_type = gro_info%ff_gromos_type
torsion_list(j)%id_type = gro_info%ff_gromos_type
@ -1343,11 +1343,11 @@ CONTAINS
imul = torsion_list(j)%torsion_kind%nmul + 1
CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul)
torsion_list(j)%torsion_kind%id_type = do_ff_charmm
torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k)
torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k)
torsion_list(j)%torsion_kind%cosphi0(imul) = chm_info%torsion_cosphi0(k)
torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k)
torsion_list(j)%torsion_kind%nmul = imul
found = .TRUE.
END IF
@ -1366,11 +1366,11 @@ CONTAINS
imul = torsion_list(j)%torsion_kind%nmul + 1
CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul)
torsion_list(j)%torsion_kind%id_type = do_ff_charmm
torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k)
torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k)
torsion_list(j)%torsion_kind%cosphi0(imul) = chm_info%torsion_cosphi0(k)
torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k)
torsion_list(j)%torsion_kind%nmul = imul
found = .TRUE.
END IF
@ -1392,11 +1392,11 @@ CONTAINS
imul = torsion_list(j)%torsion_kind%nmul + 1
CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul)
torsion_list(j)%torsion_kind%id_type = do_ff_amber
torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k)
torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k)
torsion_list(j)%torsion_kind%cosphi0(imul) = amb_info%torsion_cosphi0(k)
torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k)
torsion_list(j)%torsion_kind%nmul = imul
found = .TRUE.
END IF
@ -1415,11 +1415,11 @@ CONTAINS
imul = torsion_list(j)%torsion_kind%nmul + 1
CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul)
torsion_list(j)%torsion_kind%id_type = do_ff_amber
torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k)
torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k)
torsion_list(j)%torsion_kind%cosphi0(imul) = amb_info%torsion_cosphi0(k)
torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k)
torsion_list(j)%torsion_kind%nmul = imul
found = .TRUE.
END IF
@ -1441,11 +1441,11 @@ CONTAINS
imul = torsion_list(j)%torsion_kind%nmul + 1
CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%cosphi0,1,imul)
CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul)
torsion_list(j)%torsion_kind%id_type = inp_info%torsion_kind(k)
torsion_list(j)%torsion_kind%k(imul) = inp_info%torsion_k(k)
torsion_list(j)%torsion_kind%m(imul) = inp_info%torsion_m(k)
torsion_list(j)%torsion_kind%cosphi0(imul) = COS(inp_info%torsion_cosphi0(k))
torsion_list(j)%torsion_kind%phi0(imul) = inp_info%torsion_phi0(k)
torsion_list(j)%torsion_kind%nmul = imul
found = .TRUE.
END IF

22
src/force_fields_ext.F
View file

@ -81,7 +81,7 @@ CONTAINS
itype, iw, jatom, ncon, &
ntype, offset, stat
LOGICAL :: failure, found
REAL(KIND=dp) :: cost2, csq, sdet
REAL(KIND=dp) :: cosphi0, cost2, csq, sdet
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_parser_type), POINTER :: parser
TYPE(gromos_info_type), POINTER :: gro_info
@ -236,12 +236,13 @@ CONTAINS
CALL parser_get_object(parser,ntype,error=error)
CALL reallocate(gro_info%torsion_k,1,ntype)
CALL reallocate(gro_info%torsion_m,1,ntype)
CALL reallocate(gro_info%torsion_cosphi0,1,ntype)
CALL reallocate(gro_info%torsion_phi0,1,ntype)
DO itype=1,ntype
CALL parser_get_next_line(parser,1,error=error)
CALL parser_get_object(parser,gro_info%torsion_k(itype),error=error)
CALL parser_get_object(parser,gro_info%torsion_cosphi0(itype),error=error)
CALL parser_get_object(parser,cosphi0,error=error)
CALL parser_get_object(parser,gro_info%torsion_m(itype),error=error)
gro_info%torsion_phi0(itype) = ACOS(cosphi0)
gro_info%torsion_k(itype) = cp_unit_to_cp2k(gro_info%torsion_k(itype),"kjmol",error=error)
IF(iw>0) WRITE(iw,*) "GROMOS_FF| PUT DIHEDRALTYPE INFO HERE!!!!"
END DO
@ -523,14 +524,14 @@ CONTAINS
CALL reallocate(chm_info%torsion_d,1,ntorsion)
CALL reallocate(chm_info%torsion_k,1,ntorsion)
CALL reallocate(chm_info%torsion_m,1,ntorsion)
CALL reallocate(chm_info%torsion_cosphi0,1,ntorsion)
CALL reallocate(chm_info%torsion_phi0,1,ntorsion)
chm_info%torsion_a(ntorsion) = string
chm_info%torsion_b(ntorsion) = string2
chm_info%torsion_c(ntorsion) = string3
chm_info%torsion_d(ntorsion) = string4
CALL parser_get_object(parser,chm_info%torsion_k(ntorsion),error=error)
CALL parser_get_object(parser,chm_info%torsion_m(ntorsion),error=error)
CALL parser_get_object(parser,chm_info%torsion_cosphi0(ntorsion),error=error)
CALL parser_get_object(parser,chm_info%torsion_phi0(ntorsion),error=error)
IF(iw>0) WRITE(iw,*) " CHM TORSION ",ntorsion,&
TRIM(chm_info%torsion_a(ntorsion))," ",&
TRIM(chm_info%torsion_b(ntorsion))," ",&
@ -538,14 +539,12 @@ CONTAINS
TRIM(chm_info%torsion_d(ntorsion))," ",&
chm_info%torsion_k(ntorsion),&
chm_info%torsion_m(ntorsion),&
chm_info%torsion_cosphi0(ntorsion)
chm_info%torsion_phi0(ntorsion)
! Do some units conversion into internal atomic units
chm_info%torsion_cosphi0(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_cosphi0(ntorsion),&
chm_info%torsion_phi0(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_phi0(ntorsion),&
"deg",error=error)
chm_info%torsion_k(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_k(ntorsion),"kcalmol",&
error=error)
chm_info%torsion_cosphi0(ntorsion)=&
COS(chm_info%torsion_cosphi0(ntorsion))
CALL charmm_get_next_line(parser,1,error=error)
END DO
ELSE
@ -764,11 +763,10 @@ CONTAINS
TRIM(amb_info%torsion_d(i))," ",&
amb_info%torsion_k(i),&
amb_info%torsion_m(i),&
amb_info%torsion_cosphi0(i)
amb_info%torsion_phi0(i)
! Do some units conversion into internal atomic units
! For AMBER forcefield the cosphi0 contains the phi0 since the potential
amb_info%torsion_cosphi0(i) = cp_unit_to_cp2k(amb_info%torsion_cosphi0(i),"rad",error=error)
amb_info%torsion_phi0(i) = cp_unit_to_cp2k(amb_info%torsion_phi0(i),"rad",error=error)
amb_info%torsion_k(i) = cp_unit_to_cp2k(amb_info%torsion_k(i),"kcalmol",error=error)
END DO

14
src/force_fields_input.F
View file

@ -336,10 +336,10 @@ CPSourceFileRef,&
CALL reallocate(inp_info%torsion_d,1,ntors)
CALL reallocate(inp_info%torsion_k,1,ntors)
CALL reallocate(inp_info%torsion_m,1,ntors)
CALL reallocate(inp_info%torsion_cosphi0,1,ntors)
CALL read_torsions_section(inp_info%torsion_kind,inp_info%torsion_a, inp_info%torsion_b, inp_info%torsion_c,&
inp_info%torsion_d,inp_info%torsion_k,inp_info%torsion_cosphi0, inp_info%torsion_m,&
tmp_section, ntot, error )
CALL reallocate(inp_info%torsion_phi0,1,ntors)
CALL read_torsions_section(inp_info%torsion_kind,inp_info%torsion_a, inp_info%torsion_b,&
inp_info%torsion_c,inp_info%torsion_d,inp_info%torsion_k,inp_info%torsion_phi0,&
inp_info%torsion_m,tmp_section, ntot, error )
END IF
END IF
@ -1209,12 +1209,12 @@ CPSourceFileRef,&
!> \author teo
! *****************************************************************************
SUBROUTINE read_torsions_section(torsion_kind,torsion_a, torsion_b, torsion_c, torsion_d, torsion_k,&
torsion_cosphi0, torsion_m, section, start, error )
torsion_phi0, torsion_m, section, start, error )
INTEGER, DIMENSION(:), POINTER :: torsion_kind
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: torsion_a, torsion_b, &
torsion_c, torsion_d
REAL(KIND=dp), DIMENSION(:), POINTER :: torsion_k, torsion_cosphi0
REAL(KIND=dp), DIMENSION(:), POINTER :: torsion_k, torsion_phi0
INTEGER, DIMENSION(:), POINTER :: torsion_m
TYPE(section_vals_type), POINTER :: section
INTEGER, INTENT(IN) :: start
@ -1243,7 +1243,7 @@ CPSourceFileRef,&
CALL uppercase(torsion_c(start+isec))
CALL uppercase(torsion_d(start+isec))
CALL section_vals_val_get(section,"K",i_rep_section=isec,r_val=torsion_k(start+isec),error=error)
CALL section_vals_val_get(section,"PHI0",i_rep_section=isec,r_val=torsion_cosphi0(start+isec),error=error)
CALL section_vals_val_get(section,"PHI0",i_rep_section=isec,r_val=torsion_phi0(start+isec),error=error)
CALL section_vals_val_get(section,"M",i_rep_section=isec,i_val=torsion_m(start+isec),error=error)
END DO
END SUBROUTINE read_torsions_section

6
src/force_fields_util.F
View file

@ -976,16 +976,16 @@ CONTAINS
new_torsion_kind_set(counter)= torsion_kind_set(itorsion)
i = SIZE(torsion_kind_set(itorsion)%m)
j = SIZE(torsion_kind_set(itorsion)%k)
k = SIZE(torsion_kind_set(itorsion)%cosphi0)
k = SIZE(torsion_kind_set(itorsion)%phi0)
ALLOCATE(new_torsion_kind_set(counter)%m(i),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(new_torsion_kind_set(counter)%k(i),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(new_torsion_kind_set(counter)%cosphi0(i),STAT=stat)
ALLOCATE(new_torsion_kind_set(counter)%phi0(i),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
new_torsion_kind_set(counter)%m= torsion_kind_set(itorsion)%m
new_torsion_kind_set(counter)%k= torsion_kind_set(itorsion)%k
new_torsion_kind_set(counter)%cosphi0= torsion_kind_set(itorsion)%cosphi0
new_torsion_kind_set(counter)%phi0= torsion_kind_set(itorsion)%phi0
END IF
END DO
counter=0

8
src/input_cp2k_mm.F
View file

@ -490,10 +490,10 @@ CONTAINS
description="Define the kind of torsion potential",&
usage="KIND CHARMM",&
enum_c_vals=s2a("CHARMM","G87","G96","AMBER"),&
enum_desc=s2a("Functional Form (CHARMM|G87|G96): K * [ 1 + cos[M*PHI] * cos[PHI0]]",&
"Functional Form (CHARMM|G87|G96): K * [ 1 + cos[M*PHI] * cos[PHI0]]",&
"Functional Form (CHARMM|G87|G96): K * [ 1 + cos[M*PHI] * cos[PHI0]]",&
"Function Form AMBER : K * [ 1 + cos[M*PHI - PHI0]]"),&
enum_desc=s2a("Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]",&
"Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]",&
"Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]",&
"Functional Form (CHARMM|G87|G96|AMBER): K * [ 1 + cos[M*PHI - PHI0]]"),&
enum_i_vals=(/ do_ff_charmm,&
do_ff_g87,&
do_ff_g96,&

2
src/mc_misc.F
View file

@ -767,7 +767,7 @@ SUBROUTINE mc_make_dat_file(coordinates,natoms_tot,box_length,filename,&
TRIM(ADJUSTL(name_c)),TRIM(ADJUSTL(name_d))
WRITE(unit,20003) " M ",torsion_list(itorsion)%torsion_kind%m(imul)
WRITE(unit,20002) " K ",torsion_list(itorsion)%torsion_kind%k(imul)
WRITE(unit,20002) " PHI0 ",torsion_list(itorsion)%torsion_kind%cosphi0(imul)
WRITE(unit,20002) " PHI0 ",torsion_list(itorsion)%torsion_kind%phi0(imul)
WRITE(unit,10000) " &END TORSION"
ENDDO

95
src/mol_force.F
View file

@ -191,12 +191,13 @@ CONTAINS
!> none
!> \author DG
! *****************************************************************************
RECURSIVE SUBROUTINE force_torsions (id_type, is32, dist1, dist2, u1, u2, iw1, iw2, &
k, cosphi0, m, gt1, gt2, gt3, gt4, energy, fscalar )
SUBROUTINE force_torsions (id_type, s32, is32, ism, isn, dist1, dist2, tm,&
tn, t12, t32, k, phi0, m, gt1, gt2, gt3, gt4, energy, fscalar )
INTEGER, INTENT(IN) :: id_type
REAL(KIND=dp), INTENT(IN) :: is32, dist1, dist2
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: u1, u2
REAL(KIND=dp), INTENT(IN) :: iw1, iw2, k, cosphi0
REAL(KIND=dp), INTENT(IN) :: s32, is32, ism, isn, dist1, &
dist2
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tm, tn, t12, t32
REAL(KIND=dp), INTENT(IN) :: k, phi0
INTEGER, INTENT(IN) :: m
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4
REAL(KIND=dp), INTENT(OUT) :: energy, fscalar
@ -204,66 +205,32 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'force_torsions', &
routineP = moduleN//':'//routineN
REAL(KIND=dp) :: cosphi, lcosphi0, trig
LOGICAL :: failure
REAL(KIND=dp) :: cosphi, phi
cosphi = DOT_PRODUCT ( u1, u2 ) * iw1 * iw2
failure = .FALSE.
cosphi = DOT_PRODUCT ( tm, tn ) * ism * isn
IF(cosphi > 1.0_dp) cosphi = 1.0_dp
IF(cosphi < -1.0_dp) cosphi = -1.0_dp
phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ))
! compute energy
! Select force field
SELECT CASE (id_type)
CASE (do_ff_charmm,do_ff_g87,do_ff_g96)
SELECT CASE ( m )
CASE ( 0 )
trig = 1.0_dp
fscalar = 0.0_dp
CASE ( 1 )
trig = cosphi
fscalar = 1.0_dp
CASE ( 2 )
trig = 2.0_dp * cosphi* cosphi - 1.0_dp
fscalar = 4.0_dp * cosphi
CASE ( 3 )
trig = 4.0_dp * cosphi * cosphi * cosphi - &
3.0_dp * cosphi
fscalar = 12.0_dp * cosphi * cosphi - 3.0_dp
CASE ( 4 )
trig = 8.0_dp * cosphi * cosphi * cosphi * cosphi - &
8.0_dp * cosphi * cosphi + 1.0_dp
fscalar = 32.0_dp * cosphi * cosphi * cosphi - 16.0_dp * cosphi
CASE ( 5 )
trig = 16.0_dp * cosphi * cosphi * cosphi * cosphi * cosphi - &
20.0_dp * cosphi * cosphi * cosphi + &
5.0_dp * cosphi
fscalar = 80.0_dp * cosphi * cosphi * cosphi * cosphi - &
60.0_dp * cosphi * cosphi + 5.0_dp
CASE ( 6 )
trig = 32.0_dp * cosphi * cosphi * cosphi * cosphi * cosphi * cosphi - &
48.0_dp * cosphi * cosphi * cosphi * cosphi + &
18.0_dp * cosphi * cosphi - 1.0_dp
fscalar = 192.0_dp * cosphi * cosphi * cosphi * cosphi * cosphi - &
192.0_dp * cosphi * cosphi * cosphi + 36.0_dp * cosphi
CASE DEFAULT
CALL stop_program (routineP,"torsion parameter too large")
END SELECT
energy = k * ( 1.0_dp + cosphi0 * trig )
CASE (do_ff_charmm,do_ff_g87,do_ff_g96, do_ff_amber)
! compute energy
energy = k * ( 1.0_dp + COS(m*phi-phi0))
! compute fscalar
fscalar = -fscalar * cosphi0 * k
! compute the gradients
gt1 = ( u2 * iw2 - cosphi * u1 * iw1 ) * iw1
gt4 = ( u1 * iw1 - cosphi * u2 * iw2 ) * iw2
gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 + dist2 * is32 ** 2 * gt4
gt3 = ( - dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1
CASE (do_ff_amber)
! Use the same expression of CHARMM for AMBER (to be checked!)
lcosphi0 = COS(cosphi0)
CALL force_torsions (do_ff_charmm, is32, dist1, dist2, u1, u2, iw1, iw2, &
k, lcosphi0, m, gt1, gt2, gt3, gt4, energy, fscalar )
fscalar = k * m * SIN(m*phi-phi0)
CASE DEFAULT
CALL stop_program (routineP,"unmatched torsion kind")
END SELECT
! compute the gradients
gt1 = ( s32 * ism * ism ) * tm
gt4 = -( s32 * isn * isn ) * tn
gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4
gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1
END SUBROUTINE force_torsions
! *****************************************************************************
@ -286,14 +253,15 @@ CONTAINS
routineP = moduleN//':'//routineN
REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp
LOGICAL :: failure
REAL(KIND=dp) :: cosphi, phi
! define cosphi
failure = .FALSE.
! define cosphi
cosphi = DOT_PRODUCT ( tm, tn ) * ism * isn
IF(cosphi > 1.0_dp) cosphi = 1.0_dp
IF(cosphi > 1.0_dp) cosphi = 1.0_dp
IF(cosphi < -1.0_dp) cosphi = -1.0_dp
phi = SIGN ( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ) )
phi = SIGN( ACOS ( cosphi ), DOT_PRODUCT ( t12, tn ))
SELECT CASE (id_type)
CASE (do_ff_charmm)
@ -315,11 +283,10 @@ CONTAINS
END SELECT
! compute the gradients
gt1 = ( s32* ism * ism ) * tm
gt1 = ( s32 * ism * ism ) * tm
gt4 = -( s32 * isn * isn ) * tn
gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4
gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1
gt2 = ( dist1 * is32 ** 2 - 1.0_dp ) * gt1 - dist2 * is32 ** 2 * gt4
gt3 = ( dist2 * is32 ** 2 - 1.0_dp ) * gt4 - dist1 * is32 ** 2 * gt1
END SUBROUTINE force_imp_torsions
! *****************************************************************************

12
src/topology_amber.F
View file

@ -854,9 +854,7 @@ CPSourceFileRef,&
! ----------------------------------------------------------
! Force Fields informations related to torsions
! in amb_info%cosphi0 we store PHI0 and not its cosine, this
! is just to avoid an additional field in the torsional kind.
! Just beware of it!
! in amb_info%phi0 we store PHI0
! ----------------------------------------------------------
CALL reallocate(amb_info%torsion_a, 1,buffer_size)
CALL reallocate(amb_info%torsion_b, 1,buffer_size)
@ -864,17 +862,17 @@ CPSourceFileRef,&
CALL reallocate(amb_info%torsion_d, 1,buffer_size)
CALL reallocate(amb_info%torsion_k, 1,buffer_size)
CALL reallocate(amb_info%torsion_m, 1,buffer_size)
CALL reallocate(amb_info%torsion_cosphi0, 1,buffer_size)
CALL reallocate(amb_info%torsion_phi0, 1,buffer_size)
nsize = 0
! Torsions containing hydrogens
CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b,&
amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, &
amb_info%torsion_m, amb_info%torsion_cosphi0, particle_set, nsize,&
amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize,&
nphih, iph, jph, kph, lph, icph, pk, pn, phase, error)
! Torsions non-containing hydrogens
CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b,&
amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, &
amb_info%torsion_m, amb_info%torsion_cosphi0, particle_set, nsize,&
amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize,&
nphia, ip, jp, kp, lp, icp, pk, pn, phase, error)
! Shrink arrays size to the minimal request
CALL reallocate(amb_info%torsion_a, 1,nsize)
@ -883,7 +881,7 @@ CPSourceFileRef,&
CALL reallocate(amb_info%torsion_d, 1,nsize)
CALL reallocate(amb_info%torsion_k, 1,nsize)
CALL reallocate(amb_info%torsion_m, 1,nsize)
CALL reallocate(amb_info%torsion_cosphi0, 1,nsize)
CALL reallocate(amb_info%torsion_phi0, 1,nsize)
! ----------------------------------------------------------
! Post process of LJ parameters

4
tests/FE/regtest-1/TEST_FILES_RESET
View file

@ -14,3 +14,7 @@ Solv_alch_chng_simpl.inp
Solv_alch_chng.inp
Solv_alch_chng_dist.inp
Solv_alch_chng_simpl.inp
# numerics
Solv_alch_chng_dist.inp
Solv_alch_chng.inp
Solv_alch_chng_simpl.inp

2
tests/FE/regtest-2/TEST_FILES_RESET
View file

@ -6,3 +6,5 @@ Solv_alch_chng.inp
Solv_alch_chng.inp
# eps0
Solv_alch_chng.inp
# numerics
Solv_alch_chng.inp

3
tests/FE/regtest-3/TEST_FILES_RESET
View file

@ -10,3 +10,6 @@ Solv_alch_chng_res.inp
# eps0
Solv_alch_chng.inp
Solv_alch_chng_res.inp
# numerics
Solv_alch_chng.inp
Solv_alch_chng_res.inp

4
tests/Fist/regtest-3/TEST_FILES_RESET
View file

@ -99,3 +99,7 @@ nh3_dist.inp
nh3_dist_npt.inp
NaCl.inp
NaCl-H2O.inp
# numerics
ace_ala_nme-ambconn.inp
ace_ala_nme-amber.inp
ace_ala_nme.inp

5
tests/Fist/regtest-4/TEST_FILES_RESET
View file

@ -164,3 +164,8 @@ nh3-meta-walks_2.inp
nh3-meta-walks_2r.inp
# yet numerics
nh3-meta-walks_2r.inp
# numerics
cycbut.inp
cyhex.inp
gly_amber.inp
loop.inp

9
tests/Fist/regtest-5/TEST_FILES_RESET
View file

@ -158,3 +158,12 @@ JAC_gen.inp
# numerics after new units converter
N3dye_vib_bfgs2.inp
N3dye_vib_bfgs.inp
# numerics
arginine0.inp
arginine_crossC.inp
arginine.inp
N3dye_geoopt.inp
N3dye_vib_bfgs2.inp
N3dye_vib_bfgs.inp
N3dye_vib_restart_vec.inp
phenobenz.inp

4
tests/Fist/regtest/TEST_FILES_RESET
View file

@ -463,3 +463,7 @@ multipole_quadrupole.dbg_f_rec.inp
multipole_dipole.dbg_f_rec.inp
multipole_dip_qu.dbg_f_rec.inp
multipole_quadrupole.dbg_f_rec.inp
# numerics
multi_frag_check.inp
multi_frag.inp
multi_psf.inp

5
tests/NEB/regtest-1/TEST_FILES_RESET
View file

@ -44,3 +44,8 @@
# numerics after new units converter
2gly_CI-NEB.inp
2gly_IT-NEB.inp
# numerics
2gly_B-NEB.inp
2gly_CI-NEB.inp
2gly_EB-NEB.inp
2gly_IT-NEB.inp

7
tests/NEB/regtest-2/TEST_FILES_RESET
View file

@ -58,3 +58,10 @@
2gly_DIIS-SD-2.inp
# numerics
2gly_DIIS-SD-NEB.inp
# numerics
2gly_DIIS-DNEB.inp
2gly_DIIS-NEB.inp
2gly_DIIS-OEP-NEB.inp
2gly_DIIS-SD-2.inp
2gly_DIIS-SD-NEB.inp
2gly_DIIS-SM.inp

3
tests/NEB/regtest-3/TEST_FILES_RESET
View file

@ -13,3 +13,6 @@
2gly_IT-NEB-res.inp
# numerics
2gly_IT-NEB-res.inp
# numerics
2gly_IT-NEB.inp
2gly_IT-NEB-res.inp

3
tests/NEB/regtest-4/TEST_FILES_RESET
View file

@ -5,3 +5,6 @@
# numerics after new units converter
2gly_IT-NEB-CV.inp
2gly_IT-NEB-CV-res.inp
# numerics
2gly_IT-NEB-CV.inp
2gly_IT-NEB-CV-res.inp

3
tests/QMMM/QS/regtest-3/TEST_FILES_RESET
View file

@ -413,3 +413,6 @@ constr_hb_mm.inp
constr_hb_qm.inp
# corrected build of ppnl lists
C4H10-qmmm-gauss-0.inp
# numerics
constr_hb_mm.inp
constr_hb_qm.inp

6
tests/QMMM/SE/regtest/TEST_FILES_RESET
View file

@ -96,3 +96,9 @@ water_3_fdbg_inv.inp
water_3_walls.inp
water_3_walls_q.inp
zwitt.inp
# numerics
mol_CSVR_gen1.inp
mol_CSVR_gen2.inp
mol_CSVR_gen3.inp
mol.inp
mol_sph_cut.inp