diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_speed/input_cp2k_motion.F src/input_cp2k_motion.F
*** src_speed/input_cp2k_motion.F	2012-06-04 15:40:43.912686116 +0100
--- src/input_cp2k_motion.F	2012-06-04 15:43:18.843024903 +0100
***************
*** 132,137 ****
--- 132,166 ----
              n_keywords=0, n_subsections=1, repeats=.TRUE., required=.FALSE.,&
              error=error)
         NULLIFY(keyword)
+ 
+        CALL cp_print_key_section_create(print_key,"COLVAR_ATOMS",&
+             description="controls the output of a list of atoms. Useful for colvars.",&
+             print_level=low_print_level, common_iter_levels=1,&
+             filename="",error=error)
+        !CALL section_create(subsubsection,name="COLVAR_ATOMS",&
+        !     description=" controls the output of a list of atoms. Useful for colvars.", &
+        !     n_keywords=1, n_subsections=0, repeats=.FALSE., required=.TRUE.,&
+        !     error=error)
+        CALL keyword_create(keyword, name="LIST",&
+             description="Specifies a list of atomic indexes of the isolated ion",&
+             usage="LIST {integer}", type_of_var=integer_t, n_var=-1,&
+             required=.TRUE., repeats=.TRUE.,error=error)
+        CALL section_add_keyword(print_key,keyword,error=error)
+        CALL keyword_release(keyword,error=error)
+ 
+        CALL section_add_subsection(subsection,print_key,error=error)
+        CALL section_release(print_key,error=error)
+ 
+        !NB
+        CALL cp_print_key_section_create(print_key,"RESTRAINTS",&
+             description="controls the output of restraint values, energies forces, for colvar only right now",&
+             print_level=low_print_level, common_iter_levels=1,&
+             filename="",error=error)
+        CALL section_add_subsection(subsection,print_key,error=error)
+        CALL section_release(print_key,error=error)
+        !NB
+ 
+ 
         CALL cp_print_key_section_create(print_key,"trajectory",&
              description="controls the output of the trajectory",&
              print_level=low_print_level, common_iter_levels=1,&
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_speed/md_energies.F src/md_energies.F
*** src_speed/md_energies.F	2012-06-04 15:42:13.040053165 +0100
--- src/md_energies.F	2012-06-04 15:43:19.091965946 +0100
***************
*** 49,54 ****
--- 49,56 ----
                                               section_release,&
                                               section_type,&
                                               section_vals_get_subs_vals,&
+                                              section_vals_get,&
+                                              section_vals_val_get,&
                                               section_vals_type
    USE kinds,                           ONLY: default_string_length,&
                                               dp
***************
*** 59,70 ****
--- 61,74 ----
    USE md_environment_types,            ONLY: get_md_env,&
                                               md_environment_type,&
                                               set_md_env
+   USE memory_utilities,                ONLY: reallocate
    USE motion_utils,                    ONLY: write_simulation_cell,&
                                               write_stress_tensor,&
                                               write_trajectory
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: write_structure_data
    USE physcon,                         ONLY: angstrom,&
+                                              evolt, &
                                               femtoseconds,&
                                               kelvin
    USE qmmm_types,                      ONLY: qmmm_env_qm_type
***************
*** 77,82 ****
--- 81,87 ----
    USE timings,                         ONLY: timeset,&
                                               timestop
    USE virial_types,                    ONLY: virial_type
+   USE molecule_types_new,              ONLY: global_constraint_type
  #include "cp_common_uses.h"
  
    IMPLICIT NONE
***************
*** 307,312 ****
--- 312,325 ----
      TYPE(thermostats_type), POINTER          :: thermostats
      TYPE(virial_type), POINTER               :: virial
  
+     TYPE(global_constraint_type), POINTER    :: gci
+ 
+     TYPE(section_vals_type), POINTER         :: print_restraints_section, print_colvar_atoms_section, print_section
+     INTEGER, DIMENSION(:), POINTER           :: isolated_atoms, tmp_v, colvar_atoms
+     logical :: explicit
+     integer :: iatom, n_rep, traj_unit, stat, i_restraint, restraints_unit
+ 
+ 
      failure = .FALSE.
      NULLIFY(logger)
      logger => cp_error_get_logger(error)
***************
*** 471,476 ****
--- 484,558 ----
         CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot,&
              "FORCES", my_pos, my_act, middle_name="frc", extended_xmol_title=.TRUE., error=error)
  
+        ! Print atoms given in PRINT%COLVAR_ATOMS
+        ! get list
+        print_section => section_vals_get_subs_vals(motion_section,"PRINT",error=error)
+        print_colvar_atoms_section => section_vals_get_subs_vals(print_section,"COLVAR_ATOMS",error=error)
+        CALL section_vals_get(print_colvar_atoms_section, explicit=explicit, error=error)
+        nullify(colvar_atoms)
+        ALLOCATE(colvar_atoms(0), stat=stat)
+        CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
+        IF (explicit) THEN
+           CALL section_vals_val_get(print_colvar_atoms_section,"LIST", n_rep_val=n_rep, error=error)
+           DO i = 1, n_rep
+              CALL section_vals_val_get(print_colvar_atoms_section,"LIST", i_vals=tmp_v, i_rep_val=i, error=error)
+              CALL reallocate(colvar_atoms, 1, SIZE(colvar_atoms)+SIZE(tmp_v))
+              colvar_atoms(SIZE(colvar_atoms)-SIZE(tmp_v)+1:SIZE(colvar_atoms)) = tmp_v
+           END DO
+           ! print colvar
+           ! Get the otuput format
+ !          CALL get_output_format(root_section,"MOTION%PRINT%COLVAR_ATOMS",my_form,my_ext,error)
+           traj_unit = cp_print_key_unit_nr(logger,root_section,"MOTION%PRINT%COLVAR_ATOMS",&
+                   extension=".cv", file_position=my_pos, file_action=my_act, is_new_file=new_file, error=error) !&
+ !                  file_form=my_form, middle_name=TRIM(my_middle), is_new_file=new_file,&
+ !                  error=error)
+           IF (traj_unit > 0) THEN
+              IF (new_file) THEN
+                 WRITE (traj_unit,'("#",1X,A,3X,A,2X,A,2X,A,5X,A,5X,A,4X,A,2X,A,2X,A,2X,A,2X,A,2X,A)') &
+                      "Step Nr.","Time[fs]","Atom Nr.", &
+                      "Pos_x[A]","Pos_y[A]","Pos_z[A]",&
+                      "Vel_x[A/fs]","Vel_y[A/fs]","Vel_z[A/fs]", &
+                      "Frc_x[eV/A]","Frc_y[eV/A]","Frc_z[eV/A]"
+              END IF
+              do iatom=1,size(colvar_atoms)
+                WRITE (UNIT=traj_unit,FMT="(I8,F13.3,I8,9F13.8)") &
+                  !"CV| ", &
+                  itimes,time*femtoseconds, colvar_atoms(iatom), &
+                  particles%els(colvar_atoms(iatom))%r(1:3)*angstrom, &
+                  particles%els(colvar_atoms(iatom))%v(1:3)*angstrom/femtoseconds, &
+                  particles%els(colvar_atoms(iatom))%f(1:3)*evolt/angstrom
+              enddo
+           ENDIF
+           DEALLOCATE(colvar_atoms,stat=stat)
+           CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
+           CALL cp_print_key_finished_output(traj_unit,logger,root_section,"MOTION%PRINT%COLVAR_ATOMS",error=error)
+        END IF
+ 
+        !NB restraints
+        call cp_subsys_get(subsys=subsys, gci=gci, error=error)
+        IF (gci%nrestraint > 0) THEN
+ 	  print_restraints_section => section_vals_get_subs_vals(print_section,"RESTRAINTS",error=error)
+ 	  CALL section_vals_get(print_restraints_section, explicit=explicit, error=error)
+ 	  IF (explicit) THEN
+ 	     restraints_unit = cp_print_key_unit_nr(logger,root_section,"MOTION%PRINT%RESTRAINTS",&
+ 		     extension=".restraints", file_position=my_pos, file_action=my_act, is_new_file=new_file, error=error)
+ 	     IF (restraints_unit > 0) THEN
+ 	       IF (new_file) THEN
+ 		  WRITE (unit=restraints_unit, fmt='("#",7(2X,A))') "Step Nr.", "Time[fs]", "Restraint Nr.", &
+ 		     "Value[arb]", "Target[arb]", "Energy[eV]", "Force[eV/arb]"
+ 	       ENDIF
+ 	       do i_restraint=1, gci%ncolv%nrestraint
+ 		  write (unit=restraints_unit, fmt='(I8,X,F10.2,X,I5,4(X,F15.10))') &
+ 		     itimes, time*femtoseconds, i_restraint, &
+ 		     gci%lcolv(i_restraint)%value, gci%lcolv(i_restraint)%target_value, &
+ 		     gci%lcolv(i_restraint)%energy * evolt, gci%lcolv(i_restraint)%value_gradient * evolt
+ 	       end do
+ 	     ENDIF
+ 	     CALL cp_print_key_finished_output(restraints_unit,logger,root_section,"MOTION%PRINT%RESTRAINTS",error=error)
+ 	  ENDIF ! explicit
+        ENDIF ! gci%nrestraint > 0
+        !NB
+ 
         ! Print Simulation Cell
         CALL write_simulation_cell(cell, motion_section, itimes, time*femtoseconds, my_pos, my_act, error)
  
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_speed/molecule_types_new.F src/molecule_types_new.F
*** src_speed/molecule_types_new.F	2012-06-04 15:41:59.661572649 +0100
--- src/molecule_types_new.F	2012-06-04 15:43:19.675676241 +0100
***************
*** 44,49 ****
--- 44,50 ----
       TYPE(colvar_type),POINTER     :: colvar
       TYPE(colvar_type),POINTER     :: colvar_old
       REAL(KIND = dp)               :: lambda, sigma
+      REAL(KIND=dp)                 :: value, target_value, energy, value_gradient
    END TYPE local_colvar_constraint_type
  
  ! *****************************************************************************
diff --new-file --exclude '.#*' --exclude '*.orig' --exclude '*.rej' -C 3 src_speed/restraint.F src/restraint.F
*** src_speed/restraint.F	2012-06-04 15:42:01.458429666 +0100
--- src/restraint.F	2012-06-04 15:43:19.783276832 +0100
***************
*** 618,623 ****
--- 618,627 ----
            ind = lcolv(iconst) % colvar % i_atom (iatm)
            force(:,ind) = force(:,ind)-2.0_dp*k*tab*lcolv(iconst) % colvar % dsdr(:,iatm)
         END DO
+        lcolv(iconst)%value = tab+targ
+        lcolv(iconst)%target_value = targ
+        lcolv(iconst)%energy = k * tab**2
+        lcolv(iconst)%value_gradient = 2.0_dp * k * tab
      END DO
  
    END SUBROUTINE restraint_colv_low
