text:
import pandas as pd
import numpy as np
from . import constants as constant
from . import utils as util
import logging
from . import retrieve_data as rd
import json
import os
from . import lr_framework as lr_f

import tensorflow as tf
from tensorflow.python.keras import backend as K
from tensorflow.python.keras.models import Sequential, load_model
from tensorflow.python.keras.layers import Activation, TimeDistributed, Dense, Bidirectional, Dropout, LSTM
from tensorflow.python.keras.callbacks import Callback

INPUT_LAYERS = 1
HIDDEN_LAYER_SIZE = 8
DROPOUT = 0.2
COLUMN_NAMES = ['sin_H_t','cos_H_t','W_t','M_t-1']

RNN_COLUMNS =  [ 'W_t', 'sin_H_t', 'cos_H_t', 'M_t-1_0', 'M_t-1_1','M_t-1_unobs']
METHOD_NAME = "RNN"

SAVE_PATH = os.path.join(os.path.expanduser("~"),"Documents/PhD_Work/home_occupancy_state",METHOD_NAME,"results")
#SAVE_PATH = os.path.join(os.path.expanduser("~"),"Documents/PhD_Work/home_occupancy_state",METHOD_NAME,"results")

TEST_SPECIFIC_NAME = ''

def column_map(col_, val):
    if int(col_[-1]) == val:
        return 1
    else:
        return 0

def test_data_validity_for_train_test(test_Xarray, test_Yarray, train_Xarray, train_Yarray):
    for array_ in [test_Xarray, test_Yarray, train_Xarray, train_Yarray]:
        if array_.shape[0]==0:
            return False
    return True

def form_arrays(encode, time_step_buffer=12, time_step_ahead = 1):
    X_COLUMNS = RNN_COLUMNS
    size_ = time_step_buffer+time_step_ahead
    row = 0
    index_list = []
    X_, Y_ = np.array([]), np.array([])
    for ind_ in encode.index[time_step_buffer+time_step_ahead:]:
        df_sample = encode[ind_-pd.Timedelta('0.5H')*(size_-1):ind_].copy()

        if df_sample.shape[0] == size_:
            index_list.append(df_sample.index[-1])

            if time_step_ahead-1 > 0:
                df_sample.at[df_sample.index[-(time_step_ahead-1)]:,'M_t-1_unobs'] = 1
                df_sample.at[df_sample.index[-(time_step_ahead-1)]:,'M_t-1_0'] = 0
                df_sample.at[df_sample.index[-(time_step_ahead-1)]:,'M_t-1_1'] = 0
            if row == 0:
                X_ = df_sample.as_matrix(columns=X_COLUMNS).reshape(1,size_,len(X_COLUMNS))
                Y_ = df_sample.as_matrix(columns=['M_t']).reshape(1,size_,1)
            else:
                new_slice = df_sample.as_matrix(columns=X_COLUMNS).reshape(1,size_,len(X_COLUMNS))
                X_ = np.vstack((X_,new_slice))
                Y_ = np.vstack((Y_, df_sample.as_matrix(columns=['M_t']).reshape(1,size_,1)))

            row += 1

    return X_, Y_, index_list

def check_if_previously_completed(tstat, test_case):
    #TEST_SPECIFIC_NAME
    if os.path.isfile(os.path.join(SAVE_PATH,'training_{}_{}.json'.format(tstat,test_case))):
        return False
    else:
        return True

class rnn_test(lr_f.lr_test):
    def __init__(self, tstat_id):
        lr_f.lr_test.__init__(self, tstat_id)

    def encode_values(self):
        self.encoded = self.data.filter(items=['M_t-1','M_t','W_t'])
        self.encoded['sin_H_t'] = np.sin(2*np.pi*self.data['H_t']/48)
        self.encoded['cos_H_t'] = np.cos(2*np.pi*self.data['H_t']/48)

        for col_ in ['M_t-1_0','M_t-1_1']:
            self.encoded[col_] = self.data.filter(items=['M_t-1']).applymap(lambda x: column_map(col_, x))

        self.encoded['M_t-1_unobs'] = 0

    def generate_model(self):
        self.rnn_model = generate_uni_directional_RNN_model()

    def train_test_split(self, split_day):
        self.train, self.test = util.make_train_and_test_dataframes(self.encoded, split_day, time_step_buffer=12)

    def form_train_test_to_arrays(self):
        self.X_train, self.Y_train, self.train_index = form_arrays(self.train)
        self.X_test, self.Y_test, self.test_index = form_arrays(self.test)

    def generate_uni_directional_RNN_model(self):
        tf.reset_default_graph()
        try:
            del model
        except:
            pass
        K.clear_session()

        FEATURES = self.X_test.shape[2]
        LOOKBACK = self.X_test.shape[1]

        model = Sequential()

        for layer_number in range(INPUT_LAYERS):
            model.add(LSTM(HIDDEN_LAYER_SIZE, input_shape=(LOOKBACK, FEATURES), kernel_initializer="he_normal",
                                     return_sequences=True))
            model.add(Dropout(DROPOUT = 0.2))
        model.add(TimeDistributed(Dense(1, kernel_initializer="he_normal")))

        model.add(Activation('sigmoid'))
        model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['acc'])
        return model

    def generate_bi_directional_RNN_model(self):
        tf.reset_default_graph()
        try:
            del model
        except:
            pass

        K.clear_session()

        FEATURES = self.X_test.shape[2]
        LOOKBACK = self.X_test.shape[1]

        model = Sequential()

        model.add(Bidirectional(LSTM(HIDDEN_LAYER_SIZE, kernel_initializer="he_normal",
                                 return_sequences=True),input_shape=(LOOKBACK, FEATURES)))
        model.add(Dropout(DROPOUT = 0.2))
        model.add(TimeDistributed(Dense(1, kernel_initializer="he_normal")))
        model.add(Activation('sigmoid'))
        model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['acc'])
        return model

    def train_RNN(self, model):
        num_epoch = 20
        n_batch_size = 4
        history = model.fit(self.X_train, self.Y_train, batch_size=n_batch_size,
                            epochs=num_epoch, validation_data=(self.X_test, self.Y_test))
        self.train_accuracy =  history.history['acc']
        self.test_accuarcy = history.history['val_acc']
        return history

    def calculate_daily_error_rates(model, X_val, y_val, look_back=6):
        error_rate = []
        for ind in range(X_val.shape[0]):
            df_result = pd.DataFrame(columns = ["real","predict"])
            rowX, rowy = X_val[np.array([ind])], y_val[np.array([ind])]
            preds = model.predict_classes(rowX, verbose=0)
            df_result['real'] = rowy.flatten()
            df_result['predict'] = preds.flatten()
            df_result['error'] = df_result.apply(lambda row: map_error_fnc(row['real'], row['predict']), axis=1)
            error_rate.append(df_result.loc[6:,'error'].sum())
        return error_rate

    def populate_prediction_comparison(self, model):
        self.df_simulated_result = pd.DataFrame(index = self.test_index, columns = ["real","predict","error"])
        for ind in range(self.X_test.shape[0]):
            df_result = pd.DataFrame(columns = ["real","predict"], index=range(self.X_test.shape[1]))
            rowX, rowy = self.X_test[np.array([ind])], self.Y_test[np.array([ind])]

            preds = model.predict_classes(rowX, verbose=0)
            df_result['real'] = rowy.flatten()
            df_result['predict'] = preds.flatten()
            df_result['error'] = df_result.apply(lambda row: util.map_error_fnc(row['real'], row['predict']), axis=1)
            self.df_simulated_result.iat[ind,0] = df_result.iloc[-1,0]
            self.df_simulated_result.iat[ind,1] = df_result.iloc[-1,1]
            self.df_simulated_result.iat[ind,2] = df_result.iloc[-1,2]

    def save_test_results(self,trial):
        if not os.path.isdir(SAVE_PATH):
            os.path.mkdir(SAVE_PATH)
        self.df_simulated_result.to_csv(os.path.join(SAVE_PATH,"{}_{}_{}.csv.gz".format(METHOD_NAME,self.thermostat.tstat_id,trial)),
                                        compression='gzip')

    def save_train_data(self,trial,history):
        if not os.path.isdir(SAVE_PATH):
            os.mkdir(SAVE_PATH)
        training_eval = {}
        training_eval['train'] = history.history['acc']
        training_eval['train'] = history.history['val_acc']

        with open(os.path.join(SAVE_PATH,'training_{}_{}.json'.format(self.thermostat.tstat_id,trial)), 'w') as f:
            json.dump(training_eval, f)

    def run(self):
        logging.info("thermostat {} - main execution for RNN run".format(self.thermostat.tstat_id))
        self.encode_values()
        test_days = util.import_train_days()
        logging.info("thermostat {} - load train days".format(self.thermostat.tstat_id))

        for test_case in ['test_1','test_2','test_3','test_4']:
            logging.info("thermostat {} - beginning {}".format(self.thermostat.tstat_id, test_case))

            if check_if_previously_completed(self.thermostat.tstat_id, test_case):
                self.train_test_split(test_days[test_case])
                logging.info("thermostat {} - {} train/test".format(self.thermostat.tstat_id, test_case))
                self.form_train_test_to_arrays()

                if test_data_validity_for_train_test(self.X_test, self.Y_test, self.X_train, self.Y_train):
                    logging.info("thermostat {} - {} train/test array coded".format(self.thermostat.tstat_id, test_case))
                    logging.info("thermostat {} - {} train shape X:{} Y:{}".format(self.thermostat.tstat_id, test_case,
                                self.X_train.shape, self.Y_train.shape))
                    logging.info("thermostat {} - {} test shape X:{} Y:{}".format(self.thermostat.tstat_id, test_case,
                                self.X_test.shape, self.Y_test.shape))
                    # model initalized
                    model = self.generate_uni_directional_RNN_model()
                    logging.info("thermostat {} - {} generate model".format(self.thermostat.tstat_id, test_case))

                    history = self.train_RNN(model)
                    self.save_train_data(test_case,history)
                    logging.info("thermostat {} - {} saved training info".format(self.thermostat.tstat_id, test_case))

                    self.populate_prediction_comparison(model)
                    logging.info("thermostat {} - ran test".format(self.thermostat.tstat_id, test_case))

                    self.save_test_results(test_case)
                    logging.info("thermostat {} - saved results".format(self.thermostat.tstat_id, test_case))
                else:
                    logging.info("thermostat {} - {} invalid data aborting test".format(self.thermostat.tstat_id, test_case))
            else:
                logging.info("thermostat {} - {} Hey look we already did this one!".format(self.thermostat.tstat_id, test_case))


text:
# coding: utf-8

# In[39]:

import seren3
iout = 108
# sim = seren3.load("RT2")
sin_sim = seren3.load("BPASS_SIN")
bin_sim = seren3.load("BPASS_BIN")

sin_snap = sin_sim[iout]
bin_snap = bin_sim[iout]
cosmo = sin_snap.cosmo

import matplotlib
import matplotlib.pylab as plt
from seren3.analysis.plots import fit_scatter, fit_median
nbins=5

# In[60]:

import numpy as np
import pickle
import random
from seren3.scripts.mpi import write_fesc_hid_dict

def load_fesc(snapshot):

    db = write_fesc_hid_dict.load_db(snapshot.path, snapshot.ioutput)
    hids = db.keys()

    mvir = np.zeros(len(hids))
    fesc = np.zeros(len(hids))
    nphotons = np.zeros(len(hids))

    count=0.0
    for i in range(len(hids)):
        hid = hids[i]
        res = db[hid]

        ifesc = res["fesc"]
        if ifesc > 1. and ifesc <= 2:
            fesc[i] = random.uniform(0.9, 1.0)
            count += 1.0
        elif ifesc > 0. and ifesc <= 1.:
            fesc[i] = ifesc
        else:
    #         if (res["fesc"] > 10.):
    #             print "%e   %e" % (res["tot_mass"], res["fesc"])
            continue
        mvir[i] = res["hprops"]["mvir"]

        Nion_d_now = db[hid]["Nion_d_now"].in_units("s**-1 Msol**-1")
        star_mass = db[hid]["star_mass"].in_units("Msol")
        nphotons[i] = (Nion_d_now * star_mass).sum()

    print count/float(len(mvir))
    print len(mvir)

    ix = np.where(fesc > 0)
    fesc = fesc[ix]
    mvir = mvir[ix]
    nphotons = nphotons[ix]

    log_mvir = np.log10(mvir)
    log_fesc = np.log10(fesc)

    ix = np.where(np.logical_and(log_mvir >= 7.5, np.log10(fesc*100.) >= -1))
    # ix = np.where(log_mvir >= 7.5)
    log_mvir = log_mvir[ix]
    fesc = fesc[ix]
    nphotons = nphotons[ix]

    ix = np.where(~np.isnan(fesc))
    log_mvir = log_mvir[ix]
    fesc = fesc[ix]
    nphotons = nphotons[ix]

    print 'Loaded data for %d halos' % len(log_mvir)
    return log_mvir, fesc, None, nphotons

def load_tint_fesc(snapshot):
#     fname = "%s/pickle/ConsistentTrees/fesc_database_%05i.p" % (snapshot.path, snapshot.ioutput)
    fname = "%s/pickle/ConsistentTrees/time_int_fesc_all_halos_%05i.p" % (snapshot.path, snapshot.ioutput)
    data = pickle.load(open(fname, "rb"))

    db = write_fesc_hid_dict.load_db(snapshot.path, snapshot.ioutput)

    mvir = []
    fesc = []
    tint_fesc = []
    # hids = []
    nphotons = []

    for i in range(len(data)):
        res = data[i].result

        ifesc = res["fesc"][0]
        if ifesc > 1. and ifesc <= 10.:
            fesc.append(random.uniform(0.9, 1.0))
        elif ifesc > 0. and ifesc <= 1.:
            fesc.append(ifesc)
        else:
    #         if (res["fesc"] > 10.):
    #             print "%e   %e" % (res["tot_mass"], res["fesc"])
            continue
        mvir.append(res["Mvir"])
        tint_fesc.append(res["tint_fesc_hist"][0])
        hid = int(data[i].idx)
        # hids.append(hid)

        Nion_d_now = db[hid]["Nion_d_now"].in_units("s**-1 Msol**-1")
        star_mass = db[hid]["star_mass"].in_units("Msol")
        nphotons.append( (Nion_d_now * star_mass).sum() )

    mvir = np.array(mvir)
    fesc = np.array(fesc)
    tint_fesc = np.array(tint_fesc)
    # hids = np.array(hids)
    nphotons = np.array(nphotons)

#     print count/float(len(mvir))

    ix = np.where(fesc > 0)
    fesc = fesc[ix]
    tint_fesc = tint_fesc[ix]
    mvir = mvir[ix]
    # hids = hids[ix]
    nphotons = nphotons[ix]

    log_mvir = np.log10(mvir)
    log_fesc = np.log10(fesc)

    # ix = np.where(np.logical_and(log_mvir >= 7., log_fesc > -3))
    ix = np.where(log_mvir >= 7.5)
    log_mvir = log_mvir[ix]
    fesc = fesc[ix]
    tint_fesc = tint_fesc[ix]
    # hids = hids[ix]
    nphotons = nphotons[ix]

    ix = np.where(~np.isnan(tint_fesc))
    log_mvir = log_mvir[ix]
    fesc = fesc[ix]
    tint_fesc = tint_fesc[ix]
    # hids = hids[ix]
    nphotons = nphotons[ix]

    print 'Loaded data for %d halos' % len(log_mvir)
    return log_mvir, fesc, tint_fesc, nphotons

# data = {"BIN" : load_tint_fesc(bin_snap), "SIN" : load_tint_fesc(sin_snap)}
# # data = {"BIN" : load_fesc(bin_snap), "SIN" : load_fesc(sin_snap)}

# # In[83]:

# import matplotlib
# import matplotlib.pylab as plt
# from seren3.analysis.plots import fit_scatter, fit_median

# # plt.rcParams['figure.figsize'] = (8,8)

# fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16,6))
# ax = axs[0]

# # fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(6,6))
# # ax = axs

# nbins=5

# cols = ["r", "b"]
# binned_data = {}
# cosmo = bin_sim[iout].cosmo
# for key, c in zip(data.keys(), cols):
#     log_mvir, fesc, tint_fesc, hids = data[key]
#     # log_mvir, fesc, tint_fesc, nphotons = data[key]
#     fesc_percent = tint_fesc * 100.
#     # fesc_percent = fesc * 100.
#     bin_centres, mean, std = fit_scatter(log_mvir, fesc_percent, nbins=nbins)
#     binned_data[key] = (bin_centres, mean, std)
# #     bin_centres, median = fit_median(log_mvir, fesc_percent, nbins=nbins)
#     ax.scatter(log_mvir, fesc_percent, alpha=0.25, color=c, s=5)

# #     e = plt.errorbar(bin_centres, median, yerr=std, color=c, label='%s z=%1.2f' % (key, cosmo['z']),\
# #          fmt="o", markerfacecolor=c, mec='k', capsize=2, capthick=2, elinewidth=2, linestyle='-', linewidth=2.)

#     e = ax.errorbar(bin_centres, mean, yerr=std, color=c, label='%s z=%1.2f' % (key, cosmo['z']), fmt="o", markerfacecolor=c, mec='k', capsize=2, capthick=2, elinewidth=2, linewidth=3., linestyle="-")

# print binned_data["BIN"][1] / binned_data["SIN"][1]

# ax.legend(loc="lower right", prop={"size":16})
# ax.set_xlabel(r'log$_{10}$ M$_{\mathrm{vir}}$ [M$_{\odot}$/h]', fontsize=20)
# ax.set_ylabel(r'$\langle$f$\rangle$$_{\mathrm{esc}}$ ($\leq t_{\mathrm{H}}$) [%]', fontsize=20)
# # ax.set_ylabel(r'f$_{\mathrm{esc}}$(t) [%]', fontsize=20)
# ax.set_yscale("log")
# ax.set_ylim(1e0, 1.5e2)
# # plt.savefig("/home/ds381/fesc_bpass_bin_sin_00109.pdf", format="pdf")
# # plt.show()
# # plt.figure()

# plt.tight_layout()

# ax = axs[1]

# cols = ["r", "b"]
# binned_data = {}
# for key, c in zip(data.keys(), cols):
# #     log_mvir, fesc, tint_fesc, hids = data[key]
#     log_mvir, fesc, tint_fesc, nphotons = data[key]
#     log_nphotons_esc = np.log10(nphotons * tint_fesc)
#     bin_centres, mean, std = fit_scatter(log_mvir, log_nphotons_esc, nbins=nbins)
#     binned_data[key] = (bin_centres, mean, std)
#     # bin_centres, median = fit_median(log_mvir, log_fesc, nbins=nbins)
#     ax.scatter(log_mvir, log_nphotons_esc, alpha=0.25, color=c, s=5)

#     # e = plt.errorbar(bin_centres, median, yerr=std, color=c, label='%s z=%1.2f' % ("BPASS_BIN", cosmo['z']),\
#     #      fmt="o", markerfacecolor=c, mec='k', capsize=2, capthick=2, elinewidth=2, linestyle='-', linewidth=2.)

#     e = ax.errorbar(bin_centres, mean, yerr=std, color=c, label='%s z=%1.2f' % (key, cosmo['z']),         fmt="o", markerfacecolor=c, mec='k', capsize=2, capthick=2, elinewidth=2, linestyle='-', linewidth=2.)

# print binned_data["BIN"][1] / binned_data["SIN"][1]

# ax.legend(prop={"size":16})
# ax.set_xlabel(r'log$_{10}$ M$_{\mathrm{vir}}$ [M$_{\odot}$/h]', fontsize=20)
# ax.set_ylabel(r'log$_{10}$ $\dot{\mathrm{N}}_{\mathrm{ion}}(t)$ $\langle$f$\rangle$ $_{\mathrm{esc}}$ ($\leq t_{\mathrm{H}}$) [#/s]', fontsize=20)
# # plt.yscale("log")
# # plt.ylim(-1, 2.2)
# # plt.savefig("/home/ds381/fesc_bpass_bin_sin_00109.pdf", format="pdf")
# plt.show()

ax = plt.gca()
for iout, cols in zip([108, bin_sim.redshift(9.), bin_sim.redshift(11.)], ["r", "b", "k"]):
    # cols = ["r", "b"]
    data = {"BIN" : load_fesc(bin_sim[iout]), "SIN" : load_fesc(sin_sim[iout])}
    binned_data = {}
    cosmo = bin_sim[iout].cosmo
    c = cols
    for key, ls in zip(data.keys(), ["-", "--", ":"]):
        log_mvir, fesc, tint_fesc, hids = data[key]
        if (iout == bin_sim.redshift(9.)):
            cosmo["z"] = 9.
        if (iout == bin_sim.redshift(11)):
            cosmo["z"] = 11.
        # log_mvir, fesc, tint_fesc, nphotons = data[key]
        # fesc_percent = tint_fesc * 100.
        fesc_percent = fesc * 100.
        bin_centres, mean, std, sterr = fit_scatter(log_mvir, fesc_percent, nbins=nbins, ret_sterr=True)
        binned_data[key] = (bin_centres, mean, std)
    #     bin_centres, median = fit_median(log_mvir, fesc_percent, nbins=nbins)
        ax.scatter(log_mvir, fesc_percent, alpha=0.25, color=c, s=5, marker=".")

    #     e = plt.errorbar(bin_centres, median, yerr=std, color=c, label='%s z=%1.2f' % (key, cosmo['z']),\
    #          fmt="o", markerfacecolor=c, mec='k', capsize=2, capthick=2, elinewidth=2, linestyle='-', linewidth=2.)

        e = ax.errorbar(bin_centres, mean, yerr=std, color=c, label='%s z=%1.2f' % (key, cosmo['z']), fmt="o", markerfacecolor=c, mec='k', capsize=2, capthick=2, elinewidth=2, linestyle=ls, linewidth=3.)

    print binned_data["BIN"][1] / binned_data["SIN"][1]

ax.legend(loc="lower right", prop={"size":16})
ax.set_xlabel(r'log$_{10}$ M$_{\mathrm{vir}}$ [M$_{\odot}$/h]', fontsize=20)
# ax.set_ylabel(r'f$_{\mathrm{esc}}$ (<t$_{\mathrm{H}}$) [%]', fontsize=20)
ax.set_ylabel(r'f$_{\mathrm{esc}}$(t) [%]', fontsize=20)
ax.set_yscale("log")
ax.set_ylim(1e-1, 1.5e2)
ax.set_xlim(7.5, 10.)
# plt.savefig("/home/ds381/fesc_bpass_bin_sin_00109.pdf", format="pdf")
# plt.show()
# plt.figure()

plt.tight_layout()
plt.show()

ext: py

reason: import numpy


text:
import numpy as np

class BaseRegion(object):

    _atom = 0.3

    def __init__(self,
                 x_min,
                 x_max,
                 y_min,
                 y_max):
        self.x_min = x_min
        self.x_max = x_max
        self.y_min = y_min
        self.y_max = y_max
        '''
        ground_position x_range/atom * y_range/atom * 2 tensor
        '''
        self.ground_position_ = self.get_ground()

    def get_ground(self):
        x = np.arange(self.x_min, self.x_max + self._atom, self._atom)
        y = np.arange(self.y_min, self.y_max + self._atom, self._atom)
        ground_x, ground_y = np.meshgrid(x, y)
        ground_x = ground_x[:, :, np.newaxis]
        ground_y = ground_y[:, :, np.newaxis]
        ground_position = np.concatenate([ground_x, ground_y], axis=2)
        return ground_position

    def set_atom(self, atom):
        self._atom = atom

    def get_atom(self):
        return self._atom

max_stars_repo_name: JimyMa/UniformBSGaussianUE

max_issues_repo_name: JimyMa/UniformBSGaussianUE

max_forks_repo_name: JimyMa/UniformBSGaussianUE


text:
import logging
from copy import deepcopy
from itertools import combinations, permutations

import numpy as np
import qcelemental as qcel

from . import bend, cart, dimerfrag, oofp
from . import optparams as op
from . import stre, tors, v3d
from .exceptions import AlgError, OptError
from .v3d import are_collinear

# Functions related to freezing, fixing, determining, and
#    adding coordinates.

def connectivity_from_distances(geom, Z):
    """
    Creates a matrix (1 or 0) to describe molecular connectivity based on
    nuclear distances
    Parameters
    ----------
    geom : ndarray
        (nat, 3) cartesian geometry
    Z : list[int]
        (nat) list of atomic numbers

    Returns
    -------
    C : ndarray
        (nat, nat)

    """
    nat = geom.shape[0]
    C = np.zeros((len(geom), len(geom)), bool)
    # logger = logging.getLogger(__name__)
    for i, j in combinations(range(nat), 2):
        R = v3d.dist(geom[i], geom[j])
        Rcov = qcel.covalentradii.get(Z[i], missing=4.0) + qcel.covalentradii.get(Z[j], missing=4.0)
        # logger.debug("Trying to connect atoms " + str(i) + ' and ' + str(j) + " distance is: " +
        #            str(qcel.covalentradii.get(Z[i], missing=4.0) + qcel.covalentradii.get(Z[j], missing=4.0)))
        if R < op.Params.covalent_connect * Rcov:
            C[i, j] = C[j, i] = True

    return C

def add_intcos_from_connectivity(C, intcos, geom):
    """
    Calls add_x_FromConnectivity for each internal coordinate type
    Parameters
    ----------
    C : ndarray
        (nat, nat) matrix desribing connectivity
        see intcosMisc.connectivity_from_distances()
    intcos : list[simple.Simple]
            (nat) list of current internal coordinates (Stre, Bend, Tors)
    geom : ndarray
        (nat, 3) cartesian geometry

    """
    add_stre_from_connectivity(C, intcos)
    add_bend_from_connectivity(C, intcos, geom)
    add_tors_from_connectivity(C, intcos, geom)
    if op.Params.include_oofp or check_if_oofp_needed(C, intcos, geom):
        add_oofp_from_connectivity(C, intcos, geom)

def add_stre_from_connectivity(C, intcos):
    """
    Adds stretches from connectivity

    Parameters
    ----------
    C : ndarray
        (nat, nat)
    intcos : list[simple.Simple]
        (nat)
    Returns
    -------

    """

    # Norig = len(intcos)
    for i, j in combinations(range(len(C)), 2):
        if C[i, j]:
            s = stre.Stre(i, j)
            if s not in intcos:
                intcos.append(s)
    # return len(intcos) - Norig  # return number added

def add_h_bonds(geom, zs: list, num_atoms):
    """Add Hydrogen bonds to a fragments coordinate list
    Parameters
    ----------
    geom : np.ndarray
    zs : list
    num_atoms : int
    Returns
    -------
    list[stre.HBond]
    Notes
    -----
    Look for electronegative atoms.
    Find hydrogen atoms between covalent radii test and 2.3 Angstroms
    Check these hydrogen atoms are already bonded to an electronegative atom
    Check bond angle >= 90 degrees
    """
    logger = logging.getLogger(__name__)

    # N, O, F, P, S, Cl as proposed by Bakken and Helgaker
    electroneg_zs = [7, 8, 9, 15, 16, 17]
    # Get atom indices (within a fragment) for the electronegative atoms present and the
    # hydrogen atoms present also get
    electronegs_present = [index for index, z in enumerate(zs) if z in electroneg_zs]
    hydrogens = [index for index, i in enumerate(zs) if i == 1]

    # Some shortcuts
    min_factor = op.Params.covalent_connect
    limit = op.Params.h_bond_connect
    cov = qcel.covalentradii.get
    h_bonds = []

    for index_i, i in enumerate(electronegs_present):
        for j in hydrogens:
            if j < i:
                # do only i < j
                break
            distance = v3d.dist(geom[i], geom[j])
            covalent_thresh = min_factor * (cov(zs[index_i], missing=4.0) + cov(1, missing=4.0))
            if limit > distance > covalent_thresh:
                for k in range(num_atoms):
                    # grab k part in electronegs_present.
                    if k in electronegs_present:
                        test_angle = v3d.angle(geom[k], geom[j], geom[i])
                        if test_angle >= (np.pi / 2):
                            h_bonds.append(stre.HBond(i, j))

                            break  # Add hydrogen bond if 1 appropriate angle in connected atoms
    return h_bonds

def add_bend_from_connectivity(C, intcos, geom):
    """
    Adds Bends from connectivity

    Parameters
    ---------
    C : ndarray
        (nat, nat) unitary connectivity matrix
    intcos : list[simple.Simple]
        (nat) list of internal coordinates
    geom : ndarray
        (nat, 3) cartesian geometry
    Returns
    -------

    """

    # Norig = len(intcos)
    nat = len(geom)
    for i, j in permutations(range(nat), 2):
        if C[i, j]:
            for k in range(i + 1, nat):  # make i<k; the constructor checks too
                if C[j, k]:
                    try:
                        val = v3d.angle(geom[i], geom[j], geom[k])
                    except AlgError:
                        pass
                    else:
                        if val > op.Params.linear_bend_threshold:
                            b = bend.Bend(i, j, k, bend_type="LINEAR")
                            if b not in intcos:
                                intcos.append(b)

                            b2 = bend.Bend(i, j, k, bend_type="COMPLEMENT")
                            if b2 not in intcos:
                                intcos.append(b2)
                        else:
                            b = bend.Bend(i, j, k)
                            if b not in intcos:
                                intcos.append(b)
    # return len(intcos) - Norig

def add_tors_from_connectivity(C, intcos, geom):
    """
    Add torisions for all bonds present and determine linearity from existance of
    linear bends

    Parameters
    ----------
    C : ndarray
        (nat, nat) connectivity matrix
    intcos : list[simple.Simple]
        (nat) list of stretches, bends, etc...
    geom : ndarray
        (nat, 3) cartesian geometry
    Returns
    -------
    """

    # Norig = len(intcos)
    Natom = len(geom)

    # Find i-j-k-l where i-j-k && j-k-l are NOT collinear.
    for i, j in permutations(range(Natom), 2):
        if C[i, j]:
            for k in range(Natom):
                if C[k, j] and k != i:

                    # ensure i-j-k is not collinear; that a regular such bend exists
                    b = bend.Bend(i, j, k)
                    if b not in intcos:
                        continue

                    for l in range(i + 1, Natom):
                        if C[l, k] and l != j:

                            # ensure j-k-l is not collinear
                            b = bend.Bend(j, k, l)
                            if b not in intcos:
                                continue

                            t = tors.Tors(i, j, k, l)
                            if t not in intcos:
                                intcos.append(t)

    # Search for additional torsions around collinear segments.
    # Find collinear fragment j-m-k
    for j, m in permutations(range(Natom), 2):
        if C[j, m]:
            for k in range(j + 1, Natom):
                if C[k, m]:
                    # ignore if regular bend
                    b = bend.Bend(j, m, k)
                    if b in intcos:
                        continue

                    # Found unique, collinear j-m-k
                    # Count atoms bonded to m.
                    nbonds = sum(C[m])

                    if nbonds == 2:  # Nothing else is bonded to m

                        # look for an 'I' for I-J-[m]-k-L such that I-J-K is not collinear
                        J = j
                        i = 0
                        while i < Natom:
                            if C[i, J] and i != m:  # i!=J i!=m
                                b = bend.Bend(i, J, k, bend_type="LINEAR")
                                if b in intcos:  # i,J,k is collinear
                                    J = i
                                    i = 0
                                    continue
                                else:  # have I-J-[m]-k. Look for L.
                                    I = i
                                    K = k
                                    l = 0
                                    while l < Natom:
                                        if C[l, K] and l != m and l != j and l != i:
                                            b = bend.Bend(l, K, J, bend_type="LINEAR")
                                            if b in intcos:  # J-K-l is collinear
                                                K = l
                                                l = 0
                                                continue
                                            else:  # Have found I-J-K-L.
                                                L = l
                                                try:
                                                    val = v3d.tors(
                                                        geom[I],
                                                        geom[J],
                                                        geom[K],
                                                        geom[L],
                                                    )
                                                except AlgError:
                                                    pass
                                                else:
                                                    t = tors.Tors(I, J, K, L)
                                                    if t not in intcos:
                                                        intcos.append(t)
                                        l += 1
                            i += 1
    # return len(intcos) - Norig

# For now, let's just check for a single central atom bonded to all others
def check_if_oofp_needed(C, intcos, geom):
    logger = logging.getLogger(__name__)
    Natom = len(C)
    maxNneighbors = max([sum(C[i]) for i in range(Natom)])
    if maxNneighbors == Natom - 1:
        logger.debug("check_if_oofp_needed() is turning oofp ON")
        return True
    else:
        return False

def add_oofp_from_connectivity(C, intcos, geom):
    # Look for:  (terminal atom)-connected to-(tertiary atom)
    Nneighbors = [sum(C[i]) for i in range(len(C))]
    terminal_atoms = [i for i in range(len(Nneighbors)) if Nneighbors[i] == 1]

    # Find adjacent atoms
    vertex_atoms = []
    for T in terminal_atoms:
        vertex_atoms.append(np.where(C[T] == True)[0][0])

    for (T, V) in zip(terminal_atoms, vertex_atoms):
        if Nneighbors[V] < 3:
            pass
        # Find at least 2 other/side atoms
        side = []
        for N in np.where(C[V] == True)[0]:
            if N == T:
                pass
            else:
                side.append(N)

        if len(side) >= 2:
            try:
                val = v3d.oofp(geom[T], geom[V], geom[side[0]], geom[side[1]])
            except AlgError:
                raise OptError("Tried to add out-of-plane angle but couldn't evaluate it.")
                pass
            else:
                oneOofp = oofp.Oofp(T, V, side[0], side[1])
                if oneOofp not in intcos:
                    intcos.append(oneOofp)

    return

def add_cartesian_intcos(intcos, geom):
    """
    Add cartesian coordinates to intcos (takes place of internal coordinates)
    Parameters
    ----------
    intcos : list[simple.Simple]
        (nat) list of coordinates
    geom : ndarray
        (nat, 3) cartesian geometry
    Returns
    -------
    """

    # Norig = len(intcos)
    Natom = len(geom)

    for i in range(Natom):
        intcos.append(cart.Cart(i, "X"))
        intcos.append(cart.Cart(i, "Y"))
        intcos.append(cart.Cart(i, "Z"))

    # return len(intcos) - Norig

def linear_bend_check(o_molsys, dq):
    """
    Searches fragments to identify bends which are quasi-linear but not
    previously identified as "linear bends".
    Parameters
    ---------
    o_molsys : MOLSYS class
    dq : ndarray

    Returns
    -------
    list
        missing linear bends
    """

    logger = logging.getLogger(__name__)
    linear_bends = []

    for frag_index, frag in enumerate(o_molsys.fragments):
        for i, intco in enumerate(frag.intcos):
            if isinstance(intco, bend.Bend):
                new_val = intco.q(frag.geom) + dq[o_molsys.frag_1st_intco(frag_index) + i]
                A, B, C = intco.A, intco.B, intco.C

                # <ABC < 0.  A-C-B should be linear bends.
                if new_val < 0.0:
                    linear_bends.append(bend.Bend(A, C, B, bend_type="LINEAR"))
                    linear_bends.append(bend.Bend(A, C, B, bend_type="COMPLEMENT"))

                # <ABC~pi. Add A-B-C linear bends.
                elif new_val > op.Params.linear_bend_threshold:
                    linear_bends.append(bend.Bend(A, B, C, bend_type="LINEAR"))
                    linear_bends.append(bend.Bend(A, B, C, bend_type="COMPLEMENT"))

        missing_bends = [b for b in linear_bends if b not in frag.intcos]
        bend_report = [f"{b}, already present.\n" if b in missing_bends else f"{b}, missing.\n" for b in linear_bends]

        if missing_bends:
            logger.warning("\n\tThe following linear bends should be present:\n %s", "\t".join(bend_report))

    return missing_bends

def frozen_stre_from_input(frozen_stre_list, o_molsys):
    """
    Freezes stretch coordinate between atoms
    Parameters
    ----------
    frozen_stre_list : list
        each entry is a list of 2 atoms, indexed from 1
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for S in frozen_stre_list:
        if len(S) != 2:
            raise OptError("Num. of atoms in frozen stretch should be 2.")

        stretch = stre.Stre(S[0] - 1, S[1] - 1, constraint="frozen")
        f = check_fragment(stretch.atoms, o_molsys)
        try:
            frozen_stretch = o_molsys.fragments[f].intcos.index(stretch)
            o_molsys.fragments[f].intcos[frozen_stretch].freeze()
        except ValueError:
            logger.info("Stretch to be frozen not present, so adding it.\n")
            o_molsys.fragments[f].intcos.append(stretch)

def ranged_stre_from_input(ranged_stre_list, o_molsys):
    """
    Creates ranged stretch coordinate between atoms
    Parameters
    ----------
    ranged_stre_list : list
        each entry is a list of 2 atoms (indexed from 1) and 2 floats
    o_molsys : molsys.Molsys
        optking molecular system
    """
    for S in ranged_stre_list:
        if len(S) != 4:
            raise OptError("Num. of entries in ranged stretch should be 4.")

        stretch = stre.Stre(S[0] - 1, S[1] - 1)
        qmin = S[2] / stretch.q_show_factor
        qmax = S[3] / stretch.q_show_factor

        f = check_fragment(stretch.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(stretch)
            o_molsys.fragments[f].intcos[I].set_range(qmin, qmax)
        except ValueError:
            logger = logging.getLogger(__name__)
            logger.info("Stretch to be ranged not present, so adding it.\n")
            stretch.set_range(qmin, qmax)
            o_molsys.fragments[f].intcos.append(stretch)

def ext_force_stre_from_input(ext_force_stre_list, o_molsys):
    """
    Creates distance coordinate with external force
    Parameters
    ----------
    ext_force_stre_list : list
        each entry is a list of 2 atoms (indexed from 1), followed by a formula
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for S in ext_force_stre_list:
        if len(S) != 3:
            raise OptError("Num. of entries in ext. force stretch should be 3.")
        stretch = stre.Stre(S[0] - 1, S[1] - 1)
        f = check_fragment(stretch.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(stretch)
            o_molsys.fragments[f].intcos[I].ext_force = S[2]
        except ValueError:
            logger.info("External force stretch not present, so adding it.\n")
            stretch.ext_force = S[2]
            o_molsys.fragments[f].intcos.append(stretch)

def frozen_bend_from_input(frozen_bend_list, o_molsys):
    """
    Freezes bend coordinates
    Parameters
    ----------
    frozen_bend_list : list
        each entry is a list of 3 atoms numbers, indexed from 1
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for B in frozen_bend_list:
        if len(B) != 3:
            raise OptError("Num. of atoms in frozen bend should be 3.")

        bendFroz = bend.Bend(B[0] - 1, B[1] - 1, B[2] - 1, constraint="frozen")
        f = check_fragment(bendFroz.atoms, o_molsys)
        try:
            freezing_bend = o_molsys.fragments[f].intcos.index(bendFroz)
            o_molsys.fragments[f].intcos[freezing_bend].freeze()
        except ValueError:
            logger.info("Frozen bend not present, so adding it.\n")
            o_molsys.fragments[f].intcos.append(bendFroz)

def ranged_bend_from_input(ranged_bend_list, o_molsys):
    """
    Creates ranged bend coordinates
    Parameters
    ----------
    ranged_bend_list : list
        each entry is a list of 3 atoms, followed by 2 floats
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for B in ranged_bend_list:
        if len(B) != 5:
            raise OptError("Num. of entries in ranged bend should be 5.")

        Rbend = bend.Bend(B[0] - 1, B[1] - 1, B[2] - 1)
        qmin = B[3] / Rbend.q_show_factor
        qmax = B[4] / Rbend.q_show_factor
        f = check_fragment(Rbend.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(Rbend)
            o_molsys.fragments[f].intcos[I].set_range(qmin, qmax)
        except ValueError:
            logger.info("Frozen bend not present, so adding it.\n")
            Rbend.set_range(qmin, qmax)
            o_molsys.fragments[f].intcos.append(Rbend)

def ext_force_bend_from_input(ext_force_bend_list, o_molsys):
    """
    Creates bend coordinate with external force
    Parameters
    ----------
    ext_force_bend_list : list
        each entry is a list of 3 atoms (indexed from 1), followed by a formula
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for B in ext_force_bend_list:
        if len(B) != 4:
            raise OptError("Num. of entries in ext. force bend should be 4.")
        eBend = bend.Bend(B[0] - 1, B[1] - 1, B[2] - 1)
        f = check_fragment(eBend.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(eBend)
            o_molsys.fragments[f].intcos[I].ext_force = B[3]
        except ValueError:
            logger.info("External force bend not present, so adding it.\n")
            eBend.ext_force = B[3]
            o_molsys.fragments[f].intcos.append(eBend)

def frozen_tors_from_input(frozen_tors_list, o_molsys):
    """
    Freezes dihedral angles
    Parameters
    ---------
    frozen_tors_list : list
        each entry is list with 4 atoms (indexed from 1)
    o_molsys: molsys.Molsys
        optking molecular system
    """
    for T in frozen_tors_list:
        if len(T) != 4:
            raise OptError("Num. of atoms in frozen torsion should be 4.")

        torsAngle = tors.Tors(T[0] - 1, T[1] - 1, T[2] - 1, T[3] - 1, constraint="frozen")
        f = check_fragment(torsAngle.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(torsAngle)
            o_molsys.fragments[f].intcos[I].freeze()
        except ValueError:
            logging.info("Frozen dihedral not present, so adding it.\n")
            o_molsys.fragments[f].intcos.append(torsAngle)

def ranged_tors_from_input(ranged_tors_list, o_molsys):
    """
    Creates ranged dihedral angles from input
    Parameters
    ---------
    ranged_tors_list : list
        each entry is list with 4 atoms plus 2 floats
    o_molsys: molsys.Molsys
        optking molecular system
    """
    for T in ranged_tors_list:
        if len(T) != 6:
            raise OptError("Num. of entries in ranged dihedral should be 6.")

        torsAngle = tors.Tors(T[0] - 1, T[1] - 1, T[2] - 1, T[3] - 1)
        qmin = T[4] / torsAngle.q_show_factor
        qmax = T[5] / torsAngle.q_show_factor
        f = check_fragment(torsAngle.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(torsAngle)
            o_molsys.fragments[f].intcos[I].set_range(qmin, qmax)
        except ValueError:
            logging.info("Frozen dihedral not present, so adding it.\n")
            torsAngle.set_range(qmin, qmax)
            o_molsys.fragments[f].intcos.append(torsAngle)

def ext_force_tors_from_input(extForceTorsList, o_molsys):
    """
    Creates tors coordinate with external force
    Parameters
    ----------
    extForceTorsList : list
        each entry is a list of 4 atoms (indexed from 1), followed by a formula
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for T in extForceTorsList:
        if len(T) != 5:
            raise OptError("Num. of entries in ext. force dihedral should be 5.")
        torsAngle = tors.Tors(T[0] - 1, T[1] - 1, T[2] - 1, T[3] - 1)
        f = check_fragment(torsAngle.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(torsAngle)
            o_molsys.fragments[f].intcos[I].ext_force = T[4]
        except ValueError:
            logger.info("External force dihedral not present, so adding it.\n")
            torsAngle.ext_force = T[4]
            o_molsys.fragments[f].intcos.append(torsAngle)

def frozen_oofp_from_input(frozenOofpList, o_molsys):
    """
    Freezes out-of-plane angles
    Parameters
    ---------
    frozenOofpList : list
        each entry is list with 4 atoms (indexed from 1)
    o_molsys: molsys.Molsys
        optking molecular system
    """
    for T in frozenOofpList:
        if len(T) != 4:
            raise OptError("Num. of atoms in frozen out-of-plane should be 4.")

        oofpAngle = oofp.Oofp(T[0] - 1, T[1] - 1, T[2] - 1, T[3] - 1, constraint="frozen")
        f = check_fragment(oofpAngle.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(oofpAngle)
            o_molsys.fragments[f].intcos[I].freeze()
        except ValueError:
            logging.info("Frozen out-of-plane not present, so adding it.\n")
            o_molsys.fragments[f].intcos.append(oofpAngle)

def ranged_oofp_from_input(ranged_oofp_list, o_molsys):
    """
    Creates ranged out-of-plane angles from input
    Parameters
    ---------
    ranged_oofp_list : list
        each entry is list with 4 atoms plus 2 floats
    o_molsys: molsys.Molsys
        optking molecular system
    """
    for T in ranged_oofp_list:
        if len(T) != 6:
            raise OptError("Num. of entries in ranged out-of-plane should be 6.")

        oofpAngle = oofp.Oofp(T[0] - 1, T[1] - 1, T[2] - 1, T[3] - 1)
        qmin = T[4] / oofpAngle.q_show_factor
        qmax = T[5] / oofpAngle.q_show_factor
        f = check_fragment(oofpAngle.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(oofpAngle)
            o_molsys.fragments[f].intcos[I].set_range(qmin, qmax)
        except ValueError:
            logging.info("Frozen out-of-plane not present, so adding it.\n")
            oofpAngle.set_range(qmin, qmax)
            o_molsys.fragments[f].intcos.append(oofpAngle)

def ext_force_oofp_from_input(ext_force_oofp_list, o_molsys):
    """
    Creates out-of-plane coordinate with external force
    Parameters
    ----------
    ext_force_oofp_list : list
        each entry is a list of 4 atoms (indexed from 1), followed by a formula
    o_molsys : molsys.Molsys
        optking molecular system
    """
    logger = logging.getLogger(__name__)
    for T in ext_force_oofp_list:
        if len(T) != 5:
            raise OptError("Num. of entries in ext. force out-of-plane should be 5.")
        oofpAngle = oofp.Oofp(T[0] - 1, T[1] - 1, T[2] - 1, T[3] - 1)
        f = check_fragment(oofpAngle.atoms, o_molsys)
        try:
            I = o_molsys.fragments[f].intcos.index(oofpAngle)
            o_molsys.fragments[f].intcos[I].ext_force = T[4]
        except ValueError:
            logger.info("External force out-of-plane not present, so adding it.")
            oofpAngle.ext_force = T[4]
            o_molsys.fragments[f].intcos.append(oofpAngle)

def frozen_cart_from_input(frozen_cart_list, o_molsys):
    """
    Creates frozen cartesian coordinates from input
    Parameters
    ----------
    frozen_cart_list : list
        each entry is list with atom number, then list of 'x','y',or'z'
    o_molsys : molsys.Molsys
    """
    logger = logging.getLogger(__name__)
    for C in frozen_cart_list:
        if len(C) != 2:
            raise OptError("Num. of entries in frozen cart should be 2.")
        at = C[0] - 1
        f = o_molsys.atom2frag_index(at)  # get frag
        for xyz in C[1]:
            newCart = cart.Cart(at, xyz, constraint="frozen")
            try:
                freezing_cart = o_molsys.fragments[f].intcos.index(newCart)
                o_molsys.fragments[f].intcos[freezing_cart].freeze()
            except ValueError:
                logger.info("\tFrozen cartesian not present, so adding it.\n")
                o_molsys.fragments[f].intcos.append(newCart)

def ranged_cart_from_input(ranged_cart_list, o_molsys):
    """
    Creates ranged cartesian coordinates from input
    Parameters
    ----------
    ranged_cart_list : list
        each entry is list with atom number, then list with only 1 of
        'x','y', or 'z', then 2 floats, min and max
        so if user wants to range x, y and z coordinates, then user should
        enter three separate entries with their ranges
    o_molsys : molsys.Molsys
    """
    logger = logging.getLogger(__name__)
    for C in ranged_cart_list:
        if len(C) != 4:
            raise OptError("Num. of entries in ranged cart should be 4.")
        atom = C[0] - 1
        f = o_molsys.atom2frag_index(atom)  # get frag
        if len(C[1]) != 1:
            raise OptError("Ranged cartesian only takes 1 of x, y, or z.")
        newCart = cart.Cart(atom, C[1][0])
        qmin = C[2] / newCart.q_show_factor
        qmax = C[3] / newCart.q_show_factor
        try:
            I = o_molsys.fragments[f].intcos.index(newCart)
            o_molsys.fragments[f].intcos[I].set_range(qmin, qmax)
        except ValueError:
            logger.info("\tRanged cartesian not present, so adding it.\n")
            newCart.set_range(qmin, qmax)
            o_molsys.fragments[f].intcos.append(newCart)

def ext_force_cart_from_input(ext_force_cart_list, o_molsys):
    """
    Creates cartesian coordinate with external force
    Parameters
    ----------
    ext_force_cart_list : list
        each entry is list with atom number, then list with only 1 of
        'x','y', or 'z', then formula for external force
    o_molsys : molsys.Molsys
    """
    logger = logging.getLogger(__name__)
    for C in ext_force_cart_list:
        if len(C) != 3:
            raise OptError("Num. of entries in ext. force Cartesian should be 3.")

        atom = C[0] - 1
        if len(C[1]) != 1:
            raise OptError("External force Cartesian takes only 1 of x/y/z.")
        newCart = cart.Cart(atom, C[1][0])
        f = o_molsys.atom2frag_index(atom)
        try:
            I = o_molsys.fragments[f].intcos.index(newCart)
            o_molsys.fragments[f].intcos[I].ext_force = C[2]
        except ValueError:
            logger.info("External force Cartesian not present, so adding it.")
            newCart.ext_force = C[2]
            o_molsys.fragments[f].intcos.append(newCart)

def check_fragment(atom_list, o_molsys):
    """Check if a group of atoms are in the same fragment (or not).
    Implicitly this function also returns a ValueError for too high atom indices.
    Raise error if different, return fragment if same.
    """
    logger = logging.getLogger(__name__)
    fragList = o_molsys.atom_list2unique_frag_list(atom_list)
    if len(fragList) != 1:
        logger.error("Coordinate contains atoms in different fragments. Not currently supported.\n")
        raise OptError("Atom list contains multiple fragments.")
    return fragList[0]

def freeze_intrafrag(o_molsys):
    if o_molsys.nfragments < 2:
        raise OptError("Fragments are to be frozen, but there is only one of them")
    for F in o_molsys.fragments:
        F.freeze()

def add_constrained_intcos(o_molsys):
    # Frozen coordinates
    if op.Params.frozen_distance:
        frozen_stre_from_input(op.Params.frozen_distance, o_molsys)
    if op.Params.frozen_bend:
        frozen_bend_from_input(op.Params.frozen_bend, o_molsys)
    if op.Params.frozen_dihedral:
        frozen_tors_from_input(op.Params.frozen_dihedral, o_molsys)
    if op.Params.frozen_oofp:
        frozen_oofp_from_input(op.Params.frozen_oofp, o_molsys)
    if op.Params.frozen_cartesian:
        frozen_cart_from_input(op.Params.frozen_cartesian, o_molsys)

    # Ranged coordinates
    if op.Params.ranged_distance:
        ranged_stre_from_input(op.Params.ranged_distance, o_molsys)
    if op.Params.ranged_bend:
        ranged_bend_from_input(op.Params.ranged_bend, o_molsys)
    if op.Params.ranged_dihedral:
        ranged_tors_from_input(op.Params.ranged_dihedral, o_molsys)
    if op.Params.ranged_oofp:
        ranged_oofp_from_input(op.Params.ranged_oofp, o_molsys)
    if op.Params.ranged_cartesian:
        ranged_cart_from_input(op.Params.ranged_cartesian, o_molsys)

    # Coordinates with extra forces
    if op.Params.ext_force_distance:
        ext_force_stre_from_input(op.Params.ext_force_distance, o_molsys)
    if op.Params.ext_force_bend:
        ext_force_bend_from_input(op.Params.ext_force_bend, o_molsys)
    if op.Params.ext_force_dihedral:
        ext_force_tors_from_input(op.Params.ext_force_dihedral, o_molsys)
    if op.Params.ext_force_oofp:
        ext_force_oofp_from_input(op.Params.ext_force_oofp, o_molsys)
    if op.Params.ext_force_cartesian:
        ext_force_cart_from_input(op.Params.ext_force_cartesian, o_molsys)

    if op.Params.freeze_intrafrag:
        freeze_intrafrag(o_molsys)

def add_dimer_frag_intcos(o_molsys):
    # Look for coordinates in the following order:
    # 1. Check for 1 or list of dicts for 'interfrag_coords' in params
    # TODO: test non-equal weights
    # 2. Check 'frag_ref_atoms' keyword.  It is less flexible than 1.
    # and lower level.  We may want to remove it in the future.
    # 3. Auto-generate reference atoms.
    # TODO: move into a molsys class function?

    if op.Params.interfrag_coords != None:
        if type(op.Params.interfrag_coords) in [list, tuple]:
            for coord in op.Params.interfrag_coords:
                c = eval(coord)
                df = dimerfrag.DimerFrag.fromUserDict(c)
                df.update_reference_geometry(o_molsys.frag_geom(df.A_idx), o_molsys.frag_geom(df.B_idx))
                o_molsys.dimer_intcos.append(df)
        else:
            c = eval(op.Params.interfrag_coords)
            df = dimerfrag.DimerFrag.fromUserDict(c)
            df.update_reference_geometry(o_molsys.frag_geom(df.A_idx), o_molsys.frag_geom(df.B_idx))
            o_molsys.dimer_intcos.append(df)

    elif op.Params.frag_ref_atoms is not None:
        # User-defined ref atoms starting from 1. Decrement here.
        # Assuming that for trimers+, the same reference atoms are
        # desired for each coordinate involving that fragment.
        frag_ref_atoms = deepcopy(op.Params.frag_ref_atoms)
        for iF, F in enumerate(frag_ref_atoms):  # fragments
            for iRP, RP in enumerate(F):  # reference points
                for iAT in range(len(RP)):  # atoms
                    frag_ref_atoms[iF][iRP][iAT] -= 1

        for A, B in combinations(range(o_molsys.nfragments), r=2):
            df = dimerfrag.DimerFrag(A, frag_ref_atoms[A], B, frag_ref_atoms[B])
            df.update_reference_geometry(o_molsys.frag_geom(A), o_molsys.frag_geom(B))
            o_molsys.dimer_intcos.append(df)

    else:  # autogenerate interfragment coordinates
        # Tolerance for collinearity of ref points. Can be mad smaller, but its
        # riskier to start wth ref points the make very large angles
        col_tol = op.Params.interfrag_collinear_tol
        for A, B in combinations(range(o_molsys.nfragments), r=2):
            xyzA = o_molsys.frag_geom(A)
            xyzB = o_molsys.frag_geom(B)
            # Choose closest two atoms for 1st reference pt.
            (refA1, refB1) = o_molsys.closest_atoms_between_2_frags(A, B)
            frag_ref_atomsA = [[refA1]]
            frag_ref_atomsB = [[refB1]]
            # Find ref. pt. 2 on A.
            if not o_molsys.fragments[A].is_atom():
                for i in range(o_molsys.fragments[A].natom):
                    if i == refA1 or are_collinear(xyzA[i], xyzA[refA1], xyzB[refB1], col_tol):
                        continue
                    refA2 = i
                    frag_ref_atomsA.append([refA2])
                    break
                else:
                    raise OptError("could not find 2nd atom on fragment {:d}".format(A + 1))
            # Find ref. pt. 2 on B.
            if not o_molsys.fragments[B].is_atom():
                for i in range(o_molsys.fragments[B].natom):
                    if i == refB1 or are_collinear(xyzB[i], xyzB[refB1], xyzA[refA1], col_tol):
                        continue
                    refB2 = i
                    frag_ref_atomsB.append([refB2])
                    break
                else:
                    raise OptError("could not find 2nd atom on fragment {:d}".format(B + 1))
            # Find ref. pt. 3 on A.
            if o_molsys.fragments[A].natom > 2 and not o_molsys.fragments[A].is_linear():
                for i in range(o_molsys.fragments[A].natom):
                    if i in [refA1, refA2] or are_collinear(xyzA[i], xyzA[refA2], xyzA[refA1], col_tol):
                        continue
                    frag_ref_atomsA.append([i])
                    break
                else:
                    raise OptError("could not find 3rd atom on fragment {:d}".format(A + 1))
            # Find ref. pt. 3 on B.
            if o_molsys.fragments[B].natom > 2 and not o_molsys.fragments[B].is_linear():
                for i in range(o_molsys.fragments[B].natom):
                    if i in [refB1, refB2] or are_collinear(xyzB[i], xyzB[refB2], xyzB[refB1], col_tol):
                        continue
                    frag_ref_atomsB.append([i])
                    break
                else:
                    raise OptError("could not find 3rd atom on fragment {:d}".format(A + 1))

            df = dimerfrag.DimerFrag(A, frag_ref_atomsA, B, frag_ref_atomsB)
            df.update_reference_geometry(o_molsys.frag_geom(A), o_molsys.frag_geom(B))
            o_molsys.dimer_intcos.append(df)

        # print('end of add_dimer_frag_intcos')
        # print(o_molsys)
    return

max_stars_repo_name: psi-rking/optking

max_issues_repo_name: psi-rking/optking

max_forks_repo_name: psi-rking/optking