Weighting of RSMD colvar is modified for the case of subsystem = list (

cp2k#3818)
hfp · Dec 18, 2024 · 971aa1d · 971aa1d
1 parent e79617a
commit 971aa1d
Show file tree

Hide file tree

Showing 4 changed files with 18 additions and 17 deletions.
diff --git a/src/colvar_methods.F b/src/colvar_methods.F
@@ -799,14 +799,12 @@ RECURSIVE SUBROUTINE colvar_read(colvar, icol, colvar_section, para_env)
          CALL section_vals_get(frame_section, n_repetition=nr_frame)
 
          colvar%rmsd_param%nr_frames = nr_frame
-         ! Calculation is aborted if reference frame are less than 2
+         ! Calculation is aborted if reference frame are less than 1 and more than 2
          CPASSERT(nr_frame >= 1 .AND. nr_frame <= 2)
          CALL read_frames(frame_section, para_env, nr_frame, colvar%rmsd_param%r_ref, &
                           colvar%rmsd_param%n_atoms)
-
          ALLOCATE (colvar%rmsd_param%weights(colvar%rmsd_param%n_atoms))
          colvar%rmsd_param%weights = 0.0_dp
-
          CALL section_vals_val_get(rmsd_section, "SUBSET_TYPE", i_val=colvar%rmsd_param%subset)
          IF (colvar%rmsd_param%subset == rmsd_all) THEN
             ALLOCATE (colvar%rmsd_param%i_rmsd(colvar%rmsd_param%n_atoms))
@@ -871,6 +869,7 @@ RECURSIVE SUBROUTINE colvar_read(colvar, icol, colvar_section, para_env)
          END IF
 
          CALL section_vals_val_get(rmsd_section, "ALIGN_FRAMES", l_val=colvar%rmsd_param%align_frames)
+
       ELSE IF (my_subsection(17)) THEN
          ! Work on XYZ positions of atoms
          wrk_section => xyz_diag_section
@@ -4383,9 +4382,6 @@ SUBROUTINE mindist_colvar(colvar, cell, subsys, particles)
             nLcoord(i) = nLcoord(i) + num_n*inv_den_n
          END DO
          expnL(i) = EXP(lambda*nLcoord(i))
-!dbg
-! write(*,*) ii,nLcoord(i),expnL(i)
-!dbg
          den_Q = den_Q + expnL(i)
       END DO
       inv_den_Q = 1.0_dp/den_Q
@@ -4404,9 +4400,6 @@ SUBROUTINE mindist_colvar(colvar, cell, subsys, particles)
             rij = pbc(rpj, rpi, cell)
             r12 = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
 
-!dbg
-! write(*,*) ii,jj,rpi(1:3),rpj(1:3),rij(1:3),r12
-!dbg
             num_Q = num_Q + r12*expnL(j)
 
             sum_rij(j) = sum_rij(j) + r12
@@ -4422,10 +4415,6 @@ SUBROUTINE mindist_colvar(colvar, cell, subsys, particles)
       qfunc = num_Q*inv_den_Q
       dqfunc_dr = dqfunc_dr*inv_den_Q
       colvar%ss = qfunc
-!dbg
-!      write(*,*) ' ss ',  colvar%ss
-! stop
-!dbg
 
       DO i = 1, n_coord_from
          dqfunc_dnL(i) = lambda*expnL(i)*inv_den_Q*(sum_rij(i) - num_Q*inv_den_Q)

diff --git a/src/colvar_utils.F b/src/colvar_utils.F
@@ -570,6 +570,7 @@ SUBROUTINE post_process_colvar(colvar, particles)
 
       CHARACTER(len=3)                                   :: name_kind
       INTEGER                                            :: i, ii, j, natoms, nkinds, nr_frame, stat
+      REAL(KIND=dp)                                      :: mtot
 
       natoms = SIZE(particles)
       IF (colvar%type_id == coord_colvar_id) THEN
@@ -749,12 +750,21 @@ SUBROUTINE post_process_colvar(colvar, particles)
       END IF
 
       IF (colvar%type_id == rmsd_colvar_id) THEN
-         IF (colvar%rmsd_param%subset == rmsd_all .OR. colvar%rmsd_param%subset == rmsd_list) THEN
+         IF (colvar%rmsd_param%subset == rmsd_all) THEN
             ! weights are masses
             DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
                ii = colvar%rmsd_param%i_rmsd(i)
                colvar%rmsd_param%weights(ii) = particles(ii)%atomic_kind%mass
             END DO
+            mtot = MINVAL(colvar%rmsd_param%weights)
+            colvar%rmsd_param%weights = colvar%rmsd_param%weights/mtot
+         ELSE IF (colvar%rmsd_param%subset == rmsd_list) THEN
+            ! all weights are = 1
+            DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
+               ii = colvar%rmsd_param%i_rmsd(i)
+               colvar%rmsd_param%weights(ii) = 1.0_dp ! particles(ii)%atomic_kind%mass
+            END DO
+
          END IF
 
          IF (colvar%rmsd_param%align_frames) THEN

diff --git a/src/input_cp2k_colvar.F b/src/input_cp2k_colvar.F
@@ -1823,7 +1823,9 @@ SUBROUTINE create_colvar_rmsd_section(section)
 
       NULLIFY (keyword, subsection, subsubsection)
       CALL keyword_create(keyword, __LOCATION__, name="SUBSET_TYPE", &
-                          description="Define the subsytem used to compute the RMSD", &
+                          description="Define the subsytem used to compute the RMSD. With ALL the displacements"// &
+                          " are mass-weighted, with LIST all weights are set to 1,"// &
+                          " with WEIGHT_LIST a list of weights is expected from input.", &
                           usage="SUBSET_TYPE ALL", &
                           enum_c_vals=s2a("ALL", "LIST", "WEIGHT_LIST"), &
                           enum_i_vals=(/rmsd_all, rmsd_list, rmsd_weightlist/), &
@@ -1846,7 +1848,7 @@ SUBROUTINE create_colvar_rmsd_section(section)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="WEIGHTS", &
-                          description="Specify weights of atoms building the subset. ", &
+                          description="Specify weights of atoms building the subset. It is used only with WEIGHT_LIST ", &
                           usage="weightS {real} {real} ..", repeats=.TRUE., &
                           n_var=-1, type_of_var=real_t)
       CALL section_add_keyword(section, keyword)

diff --git a/tests/QS/regtest-gpw-2-3/TEST_FILES b/tests/QS/regtest-gpw-2-3/TEST_FILES
@@ -19,5 +19,5 @@ H2O_meta_langevin.inp                                  1      2e-07
 Au13ico_mtd.inp                                        1      5e-13            -433.32355790139007
 # rmsd AB
 Au12_rmsd_AB_mtd.inp                                   1      3e-11            -397.79481882168199
-Au12_rmsd_A_mtd.inp                                    1      2e-11            -397.77606915975093
+Au12_rmsd_A_mtd.inp                                    1      2e-11            -397.77606999996829
 #EOF